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FOREWORD 


The  results  of  high-speed  tensile-impact  experiments  performed 
in  the  Pioneering  Research  Laboratory  have  in  the  past  been  analyzed 
by  graphical  methods  based  on  the  method  of  characteristics.  The 
present  report  is  the  first  of  a  series  of  three  in  which  a  new 
numerical  method  suitable  for  use  in  high-speed  computers  is  des¬ 
cribed. 

The  division  into  three  reports  is  based  on  the  successive 
stages  in  the  mathematical  analysis  as  carried  out  by  Professor 
Crout.  In  the  present  report  the  basic  mathematical  technique  is 
presented  and  applied  to  a  linear  material,  so  that  the  method  can 
be  tested  in  a  situation  where  the  solution  is  already  known.  In 
the  second  report  the  method  will  be  extended  to  non-linear  mater¬ 
ials,  and  in  the  third  report  the  case  of  non-linear  materials 
exhibiting  time-dependent  response  (such  as  creep)  will  be  treated. 

The  project  reference  number  for  this  work  is  lT06ll02BllA-07. 
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ABSTRACT 


* 


An  integral-equation,  modified  successive-substitution  method 
of  solving  the  equation  of  motion  of  a  linear  elastic  string  is 
derived.  The  method  is  an  extension  to  two  independent  variables 
of  a  method  previously  employed  with  only  a  single  independent 
variable.  The  solution  is  proved  to  be  correct  in  the  linear 
case  and  is  expected  to  be  valid  also  for  nonlinear  cases  and 
for  situations  in  which  the  response  of  the  material  is  time 
dependent.  Methods  of  numerical  approximation  based  on  a  lattice 
of  points  in  the  time-position  plane  are  developed,  suitable  for 
digital-computer  programming.  Typical  results  of  computer  runs 
on  linear  materials  are  presented. 
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RESPONSE  OF  POLYMERS  TO  TENSILE  IMPACT. 

1.  An  Integral -Equation,  Successive -Substitution  Solution  for  Use  in 
Digital  Computers,  and  Its  Application  to  a  Linear  Elastic  Model. 


1.  History  and  Scope 

The  results  now  being  reported  have  come  from  the  combined  efforts 
of  several  individuals  and  groups.  Professor  Prescott  D.  Crout, 
Mathematics  Department,  Massachusetts  Institute  of  Technology,  is  res¬ 
ponsible  for  the  mathematical  developments.  Mr.  Malcolm  N.  Pilsvorth, 
Jr.,  of  the  Thermodynamics  Group,  U.  S.  Army  Natick  Laboratories,  is 
responsible  for  the  experimental  data.  The  Computer  Branch  of  the 
Natick  Laboratories  did  the  detailed  programming  and  carried  out  the 
machine  computations.  Much  of  the  planning  of  the  work  was  done 
Jointly  by  Professor  Crout,  Mr.  Pilsworth,  and  Dr.  Harold  J.  Hbge  of 
the  Thermodynamics  Group. 

The  project  has  been  concerned  with  the  response  of  polymers  such 
as  nylon  to  tensile  impact  under  conditions  in  which  both  the  strain 
and  the  rate  of  strain  are  large,  and  conventional  methods  of  analy¬ 
sis  are  inadequate. 

The  present  report  describes  the  situation  when  the  work  was 
begun,  explains  the  mathematical  method  that  was  developed,  and  shows 
the  results  of  applying  the  new  method  to  strain  propagation  in  a 
linear  material  (one  having  a  stress-strain  curve  of  constant  slope) 
with  no  creep  or  other  time  dependence.  The  new  method  is  an  itera¬ 
tive  procedure  in  which  the  values  are  calculated  at  a  network  or 
lattice  of  points  in  the  position-time  plane.  The  initial  conditions 
and-  the  impact  velocity  give  the  values  of  desired  quantities  at 
certain  of  these  lattice  points;  the  propagation  of  the  stress  and 
strain  into  the  lattice  is  given  by  the  iterative  procedure. 

Since  the  solution  of  the  linear  problem  is  known,  the  present 
report  is  essentially  a  demonstration  of  the  validity  of  the  new 
method.  In  subsequent  reports  the  method  is  applied  to  problems  of 
greater  complexity,  including  some  for  which  satisfactory  solutions 
have  not  been  previously  available. 
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2.  Introduction 


* 


A  satisfactory  theory  of  the  response  of  materials  to  loading 
should  permit  us  to  calculate  the  stresses  and  strains  in  a  material,  „ 

for  any  rate  of  loading,  from  a  few  fixed  properties  or  parameters  of 
the  material.  There  is  some  latitude  in  the  choice  of  these  fixed 
properties  or  parameters  and  those  appropriate  for  small  strains  are 
not  necessarily  the  most  convenient  when  large  strains  are  to  be 
dealt  with. 

Linear  viscoelastic  theory  is  based  on  the  concepts  of  Hookean 
"springs"  and  Newtonian  "dashpots".  It  leads  to  a  complex  modulus 
of  elasticity  in  which  the  real  part  is  called  the  storage  modulus 
and  the  imaginary  part  the  loss  modulus.  It  accounts  for  both  elastic 
and  inelastic  deformations  including  creep,  relaxation,  and  recovery. 

It  holds  for  slow  and  rapid  deformations,  including  continuous  oscilla¬ 
tions.  It  fails  for  large  strains,  however,  whether  these  are  produced 
slowly  or  rapidly,  because  the  properties  of  a  material  at  large  strains 
are  no  longer  linear. 

Large  strains  produced  at  slow  or  moderate  rates  may  be  described 
by  stress-strain  curves.  A  family  of  curves,  each  one  appropriate  to 
a  different  rate  of  straining,  will  normally  be  required.  Creep  or 
relaxation  curves  may  also  be  useful;  the  principle  of  time-temperature 
equivalence  is  sometimes  useful  when  curves  taken  at  slow  rates  of  a 

strain  are  to  be  used  at  higher  rates  of  strain. 

When  strains  are  produced  rapidly,  inertial  effects  become  im¬ 
portant.  The  strains  can  no  longer  be  considered  as  uniformly  dis¬ 
tributed  throughout  a  sample;  they  propagate  as  pulses  (wave-fronts) 
that  are  subject  to  reflection  and  attenuation.  When  inertial  effects 
begin  to  dominate,  it  is  natural  to  try  a  treatment  that  emphasizes 
these  effects.  One  way  to  do  this  is  by  formally  as suming  that  stress 
is  a  function  of  strain  but  is  independent  of  time.  With  this  simpli¬ 
fication  the  propagation  of  stress  and  strain  within  a  material  may  be 
calculated  by  the  method  of  characteristics,  as  has  been  shown  by 
Von  Karman  [1],  Taylor  [2],  and  Rachmatulin  [5].  This  method  of 
attack  was  first  applied  to  metals  and  has  achieved  some  success. 

Since  time-dependent  effects  are  much  more  obvious  in  polymers  than 

in  metals,  we  should  perhaps  expect  the  rate- independent  theory  to  be  ^ 

less  satisfactory  for  polymers.  However,  when  the  strain-rate  is  high 
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it  is  quite  conceivable  that  the  stress  would  be  almost  entirely  deter 
mined  by  the  strain  and  would  depend  only  slightly  on  the  time  elapsed 
after  the  strain  had  been  created. 

The  rate- independent  theory  as  developed  by  Von  Karman  has  been 
applied  to  polymers  in  these  Laboratories  and  elsewhere,  but  it  is  not 
very  satisfactory.  Pilsworth  and  Hoge  [4],  and  Pilsworth  [5]  have 
reported  measurements  in  which  rubber  and  nylon  were  subjected  to 
tensile  impact.  In  a  typical  experiment  on  nylon  a  sample  50  cm  long 
was  impacted  at  90  m  sec  1.  Taking  2400  m  sec  1  as  the  average 
velocity  of  propagation  of  strain  in  nylon,  the  strain  pulse  produced 
by  the  impact  had  a  magnitude  of  about  90/2400  =  3-75  percent.  This 
is  1/4  of  the  breaking  strain,  which  is  about  15  percent  for  nylon. 

The  strain  pulse  traversed  the  yarn  several  times;  during  the  third  re 
flection  the  yarn  broke.  Method-of- characteristics  calculations  were 
made  for  the  experiment,  using  several  different  stress-strain  curves. 
The  results  showed  that  any  one  of  several  stress-strain  curves  would 
predict  the  observed  strain  pattern  roughly  but  there  was  no  reason¬ 
able  stress- strain  curve  that  would  explain  the  observations 
accurately. 

In  a  typical  experiment  on  rubber,  Pilsworth  and  Hoge  impacted 
a  50  cm  length  at  42  m  sec  1.  Taking  50  m  sec  1  as  the  average 
velocity  of  propagation  of  strain  in  rubber,  the  strain  pulse  pro¬ 
duced  by  the  impact  had  a  magnitude  of  about  50/42  =  119  percent. 

This  is  a  little  less  than  1/4  the  breaking  strain,  which  for  the 
type  of  rubber  used  is  about  500  percent.  The  strain  pulse  traversed 
the  sample  several  times;  breakage  occurred  during  the  third  reflec¬ 
tion  as  it  had  for  nylon.  Method-of-characteri sties  calculations  were 
made,  again  using  several  different  stress-strain  curves.  Again  the 
conclusion  was  reached  that  there  was  no  reasonable  stress-strain 
curve  that  would  give  an  accurate  representation  of  the  observed 
strain  pattern  on  the  basis  of  time -independent  theory.  The  strain 
patterns  obtained  for  rubber  as  described  above  were  modified  by 
making  a  correction  for  creep.  This  improved  the  agreement  between 
observation  and  calculation  but  the  modification  was  tedious  and  not 
very  elegant. 

In  both  of  the  experiments  cited  above,  the  strain  pulse  was 
roughly  1/4  of  the  breaking  strain,  although  the  pulse  in  the  rubber 
was  32  times  as  large  as  the  pulse  in  the  nylon.  Both  of  these  pulses 
we  characterize  as  large,  because  each  is  an  appreciable  fraction  of 
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the  breaking  strain.  The  ,tbest"  stress-strain  curve  for  use  in  each 
of  the  experiments  described  above  was  a  derived  curve,  not  a  directly 
measured  one.  It  was  selected  simply  because  it  gave  better  results 
than  other  curves  in  the  calculation  of  strain  patterns.  The  curves 
found  in  this  way  are  steeper  and  have  a  greater  curvature  than 
stress-strain  curves  measured  conventionally  at  low  rates  of  strain. 
The  evidence  is  fairly  strong  that  time -dependence  should  not  be 
ignored,  and  that  when  it  is  ignored  the  stress-strain  curve  is  dis¬ 
torted  to  partially  compensate  for  it. 

The  work  of  Smith  and  co-workers  also  establishes  the  importance 
of  time-dependence  at  high  rates  of  strain.  For  example.  Smith  and 
Fenstermaker  [6]  studied  the  longitudinal  strain  in  rubber  subjected 
to  transverse  impact.  They  found  significant  creep  at  the  point  of 
impact  during  the  first  millisecond.  During  the  next  7  milliseconds 
creep  effects  were  small,  but  the  creep  occurring  after  8  milliseconds 
could  not  be  ignored. 

Evidence  such  as  that  given  above  brought  us  to  the  conclusion 
that  a  time -independent  theory  was  inadequate  for  polymers  and  that 
some  method  should  be  found  to  include  time-dependence  from  the 
start.  Such  a  method  is  presented  in  a  series  of  three  reports,  of 
which  this  is  the  first. 

From  a  mathematical  standpoint,  the  basic  idea  consists  of 
transforming  the  relevant  partial  differential  equations  into 
integral  equations,  and  then  solving  these  by  a  modified  method  of 
successive  substitutions.  The  resulting  method  does  not  require 
linearity,  and  is  applicable  to  systems  containing  nonlinear  and 
time-dependent  materials.  In  the  present  report,  however,  it  is 
applied  only  to  the  case  where  a  stress-strain  curve  exists  and  is 
linear,  the  purpose  being  to  develop  the  integral  equation  method, 
and  check  the  accuracy  of  the  results  which  it  gives  by  comparing 
them  with  the  corresponding  exact  results  obtained  by  other  methods. 
This  report  also  contains  some  preliminary  work  which  is  relevant 
to  time  dependent  systems. 

Since  a  clear  picture  of  the  overall  problem  must  await  the 
appearance  of  the  second  and  third  reports  of  this  series,  it  should 
be  kept  in  mind  that  our  primary  objective  is  not  to  devise  a  new 
method  for  solving  linear  problems  or  for  solving  those  nonlinear 
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problems  which  can  be  solved  by  other  methods,  but  to  devise  a  method 
capable  of  solving  nonlinear  problems  involving  time -dependent  mater¬ 
ials.  Aside  from  this,  however,  the  integral  equation  method  developed 
is  new,  should  be  a  useful  tool  in  solving  other  problems,  and  hence 
should  itself  be  of  interest. 

The  integral  equation  approach  was  chosen  because  integral 
equation  methods  have  given  good  results  in  the  past  when  used  to 
solve  problems  involving  one  or  more  ordinary  differential  equations. 

In  those  problems  integral  equations  were  used  in  preference  to  dif¬ 
ference  equations  because  of  the  greater  accuracy  of  numerical  in¬ 
tegration  as  compared  with  numerical  differentiation.  Tt  is  also 
recognized  that  methods  based  on  characteristics  have  been  applied  to 
somewhat  similar  problems  [7],  and  could  conceivably  lead  to  useful 
solutions  of  the  present  problem. 

3.  Statement  of  Problem.  Preliminary  Considerations. 

The  present  investigation  pertains  to  the  material  properties  of 
nylon,  rubber,  and,  perhaps,  other  substances.  The  first  problem  which 
we  shall  consider  is  that  of  determining  the  motion  of  a  string  of 
length  L  and  cross-section  area  A.  In  the  experimental  work,  a  typi¬ 
cal  value  of  L  was  50  cm.  In  the  case  of  rubber,  the  cross  section 
was  square,  about  1/8  inch  x  1/8  inch;  in  the  case  of  nylon,  it  was 
circular,  with  a  diameter  of  about  .05  cm  (1050  denier).  This  string 
is  fixed  at  the  left  end,  and  is  initially  at  rest.  At  time  t  =  0 
the  right  end  is  given  a  velocity  vQ  toward  the  right,  as  indicated 
in  Fig.  1.  For  the  present  vQ  will  be  held  constant.  In  order  to 


Fig.  1  String  subjected  to  tensile  impact. 
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calculate  the  subsequent  motion  of  the  string,  it  is  necessary  that 
the  elastic  properties  of  the  material  be  specified  completely.  If 
they  are  not  known  they  can  be  assumed,  and  the  validity  of  the 
assumptions  checked  by  comparing  the  resulting  calculated  motion  with 
that  determined  experimentally. 

Recalling  the  distinction  between  plane  stress  and  plane  strain 
in  elasticity  theory,  we  note  that  the  present  problem  involves 
one-dimensional  stress,  but  not  one-dimensional  strain.  A  slab  of 
material  constrained  to  have  a  fixed  cross-section  would,  on  the 
contrary,  if  given  a  constant  normal  stress  over  its  surfaces  have 
a  state  of  one-dimensional  strain,  but  not  one-dimensional  stress. 

Denoting  the  normal  stress  along  the  axis  of  the  string  by  a, 
and  the  normal  strain  along  this  axis  by  e,  we  have  as  a  first 
approximation  the  relation 

a  =  Ee,  (1) 

given  by  Hooke's  law,  E  being  the  constant  "modulus  of  elasticity". 
The  strains  perpendicular  to  the  axis  of  the  string  may  be  large, 
but  it  will  not  be  necessary  to  treat  them  explicitly.  They  play  a 
role  in  determining  the  value  of  E.  Although  (l)  is  adequate  for 
many  purposes  we  wish  to  consider  and  to  devise  mathematical 
methods  for  using  more  accurate  stress-strain  relationships.  Let 
us  consider  a  few  possible  improvements  in  (l). 

First,  any  viscous  opposition  to  the  establishment  of  the 
strain  system  would  give  rise  to  a  stress  contribution  which  is 
proportional  to  the  time  rate  of  change  of  strain,  thus 

a  =  ae  +  b  (2) 

dt 

where  a  and  b  are  constants.  This  relation  would  correspond  to  a 
spring  paralleled  by  a  viscous  damper,  the  latter  being  indicated 
by  a  dashpot.  Fig.  2. 

Second,  the  coefficient  E  in  (l)  may  not  be  constant,  but  may 
depend  upon  e,  thus 


a  =  E(e)  e . 


(3) 
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Fig.  2  Viscoelastic  model  consisting  of  a  spring  paralleled  by  a 
dashpot. 


y'-§#P—  f* 

-^rrrt — r 

Fig.  5  Model  in  which  the  modulus  E  increases  with  strain. 
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It  follows  that  the  stress-strain  relationship  is  no  longer  linear. 
If  the  stress- strain  curve  is  concave  upward  it  may  be  convenient 
to  consider  this  relation  to  correspond  to  a  set  of  springs  which 
are  in  parallel,  and  are  engaged  successively  as  e  is  increased. 
Fig.  3. 


Third,  we  may  have  a  situation  wherein  an  applied  stress  gives 
rise  to  a  strain  which  is  composed  of  two  contributions:  one  of  the 
kind  indicated  by  (l),  and  the  other  of  the  kind  indicated  by  (2). 

We  thus  have 

€  =  £i  +  e2, 

a  —  ae^,  (4) 

de2  • 

o  —  be2  +  o  ^2  =  ^  when  t  —  Oj 

which  corresponds  to  two  springs  in  series,  one  being  paralleled  by 
a  viscous  damper.  Fig.  4.  In  this  case  a  suddenly  applied  stress 
gives  rise  to  an  immediate  proportional  strain  plus  a  subsequent 
creep. 

Fourth,  it  is  more  than  likely  that  the  representation  of  the 
stress-strain  relationship  can  be  improved  by  using  more  compli¬ 
cated  systems,  whose  equivalent  mechanical  network  consists  of 
some  series-parallel  or  even  more  complicated  arrangement  of  springs 
and  viscous  dampers.  Fig.  5#  for  example. 

Corresponding  to  Fig.  5  we  have 
a  =  ax  +  a2, 

€  =  +  €2  =  €3  +  €4, 

ai  =  a€l> 

de2 

a,  =  b€2  +  c  - — ,  €2  =  0  when  t  =  0, 

1  dt 

a2  =  ee3> 
de. 

a2  =  f  e4=0  when  t  =  0, 

where  a,  b,  ...,  f  are  constants.  Finally  we  note  that  the 
mechanical  circuit  elements  need  not  be  lumped  -  they  may  be 


% 
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Fig.  4  Three-element  viscoelastic  model. 


Fig.  5  Multi-element  viscoelastic  model. 
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distributed,  as  in  the  mechanical  analog  of  an  electric  cable. 


v 

The  difficulties  which  arise  in  applying  a  nonlinear  relation¬ 
ship  such  as  (3)  will  be  considered  later.  At  present  it  will  be 
mentioned  only  that  a  nonlinear  problem  can  often  be  treated  as  a 
sequence  of  linear  problems. 

We  are  now  in  a  position  to  foresee  another  kind  of  difficulty, 
which  has  nothing  to  do  with  nonlinearity.  This  arises  from  the 
fact  that  the  stress- strain  relationship  involves  not  only  the 
present  values  of  stress  and  strain,  but  also  the  past  history  of 
the  material.  More  explicitly  we  see  that  whereas  in  the  case  of 
(l)  or  (3)  the  present  value  of  the  strain  determines  the  stress, 
and  vice  versa,  this  is  not  true  in  the  other  cases  which  we  have 
just  considered.  For  example,  in  the  case  of  (2)  if  a  constant 
stress  a  be  applied,  the  strain  changes  from  its  initial  value 
zero  to  its  final  value,  but  not  abruptly.  The  difference  between 
the  two  values  dies  out  exponentially.  Conversely  the  present 
value  of  e  is  not  sufficient  to  determine  0;  de/dt  is  also 
required.  Similar  remarks  apply  in  the  case  of  (4),  (5),  and 
more  complice  ted  cases.  This  situation  is  sufficient  to  explain, 
and  to  allow  the  inclusion  of  the  phenomena  of  hysteresis  and  creep. 

Stated  briefly,  an  elementary  piece  of  material  behaves  not  like  a 
single  element  in  a  mechanical  network,  but  like  a  mechanical  system. 

In  our  problem,  a  and  6  are  functions  of  time.  It  is  of 
interest  to  take  the  Laplace  transforms  of  the  equations  for  the 
models  we  have  been  considering,  since  all  show  an  interesting 
Hooke's  law  formalism  after  transformation.  We  shall  suppose  that 
o  and  e  are  both  identically  zero  before  the  disturbance.  Denoting 
the  Laplace  transform  of  any  quantity  by  a  bar  over  that  quantity 
we  obtain  from  (l)  by  taking  the  Laplace  transform 

(6) 

+  bpi  , 

+  bp)e  .  (7)  ^ 


a  =  El 

Similarly  (2)  gives 
0  =  ae 
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In  the  case  of  (4)  we  obtain 

* 

1  »  ii  +  i2  , 

o  =  all  > 

a  =  be2  +  cpe2  j 

whence 

2  _  0  .  o 

T  'b~i~cp  ' 

or 

5  =  a.lLj-  CP?  g  . 

a  +  b  +  cp 

Finally,  (5)  becomes 

a  =  oi  +  a2  , 
i  *  ii  +  ig  =  i3  +  i4  , 

ffi  »  , 

Oi  =  (b  +  cp)  €g  , 

02  =  ei3  , 

02  =  fp€4  . 


(8) 


It  follows  that 


—  vj,  U  - 

e  =  -1  + - 1 — 

a  b  +  cp 


Og_ 

e 


Jg 

fp 


-  _  a  (b  +  cp)  - 
ai  a  +  b  +  cp  € 


n  -  efP  2 
2  -  e  +  fp 


whence 


5  =  a_jb-+  cp)  +  - 

a+b+cp  e+fp 


€ 


(9) 


It  can  now  be  seen  that  (6),  (7),  (8),  and  (9)  are  all  of  the 
form  of  Hooke's  law;  but  now  the  proportionality  factor  between 
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a  and  e  is  a  function  of  p.  More  generally  for  any  linear  system 
let  E(p)  denote  the  Laplace  transform  of  the  stress  which  results 
from  a  unit  instantaneous  strain  impulse  applied  at  t  =  0,  then 
for  any  applied  stress-time  function 


a  =  E(p)i  ; 


(10) 


(6),  (7),  (8),  and  (9)  are  special  cases  of  this  more  general  re¬ 
lationship.  The  Laplace  transform  approach  to  the  problem  will  not 
be  developed  further. 

If  instead  of  working  with  Laplace  transforms  we  work  directly 
with  a  and  6,  we  should  have  to  deal  with  the  stress-strain  relation¬ 
ship 


(11) 


where  E(t)  is  the  stress  due  to  a  unit  strain  impulse,  and  is  the 
inverse  transform  of  E(p).  Equation  (11 )  is  evidently  more  compli¬ 
cated  than  (10 ),  and  consists  of  a  superposition  integral. 

4.  Partial  Differential  Equation  Approach. 

Our  general  procedure  will  be  to  examine  the  available  mathe¬ 
matical  methods,  and,  starting  with  simpler  problems,  try  to  solve 
problems  of  increasing  difficulty.  To  have  a  definite  starting 
point,  let  us  choose  the  simplest  stress- strain  relationship  (l), 
namely 


c  =  Ee 


where  E  is  constant,  and  apply  it  to  the  system  shown  in  Fig.  1. 
In  this  system  let 

L  m  length  of  material  unstretched,  A 

A  =  cross  section  of  material  unstretched,  / 

P  =  density  of  material  unstretched,  / 

vQ  =  constant  velocity  of  right  end  of  material,  \ 

x  =  distance  from  left  end,  unstretched,  / 

s  =  displacement  of  the  particle  of  material  which  was  1 
initially  (in  the  unstretched  state)  at  distance  x  \ 

from  the  left  end,  s  is  measured  positive  to  the  1 

right.  J 
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In  setting  up  the  partial  differential  equation  we  consider  the 
translator  equilibrium  of  an  infinitesimal  interval  dx  of  the 
unstretched  piece  of  material,  Fig.  6.  Noting  (1)  and  the  fact 


*  X+dX 


Fig.  6  Infinitesimal  element  of  impacted  string. 


that 


(13) 


s  being  a  function  of  x  and  the  time  t,  we  see  that  the  force  acting 
to  the  left  on  the  left  face  of  the  element  is  AE  .  That  acting 
to  the  right  on  the  right  face  of  the  element  is  g&en  by  a  similar 
expression,  the  only  difference  being  that  ~  must  be  evaluated  at 
(x  +  dx)  instead  of  at  x.  Since 


=  !f  +!#dx 


(1U) 


good  through  first  order  terms  in  dx,  it  follows  that  the  net  trans¬ 
lator  force  on  the  element  that  arises  because  of  the  stretching  is 

“  &  ♦  g**)  -  «  gf  -  g-|  <*•  (15) 
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For  dynamic  equilibrium  this  force  must  equal  the  mass  of  the 
element  times  its  acceleration,  or 


PA  dx 


a2s 

dt2 


(16) 


Equating  (15)  and  (16)  and  canceling  Adx,  we  obtain  finally 

d2s  _  1_  a2 s 

5x2  c£  ^t2 


where 


(17) 


(18) 


(17)  is  the  one-dimensional  wave  equation,  the  wave  velocity  being 
c.  The  general  solution  is 


s  =  f  (x  -  ct)  +  g  (x  +  ct)  , 


(19) 


which  consists  of  two  waves  of  arbitrary  shape  traveling  with  velocity 
c,  one  toward  the  right  and  one  toward  the  left  [8] . 

In  order  to  apply  (19)  to  the  present  problem  we  shall  consider 
x  to  extend  from  •«  to  ®,  and  shall  determine  f  and  g  so  that 


s=0,  |-OifO<x<L,  t«0, 

s  =  0  if  x  =  0,  and  =  v  if  x  =  L  for  any  t 

dt  0 


} 


(20) 


Equation  (19)  then  gives  the  desired  solution  for  0  x  <  L.  The  first 
two  relations  in  (20)  give 


0  =  f  (x ) 
0  =  -cf 


)  +  g(x)  ,  1 

7°  < 

'  (x)  +  eg'  (x)  , / 


x  <  L, 


(21) 


where  the  primes  indicate  differentiation.  Canceling  c  and  inte¬ 
grating  the  second  of  these  relations,  (21)  becomes 


0  =  f  (x )  +  g  (x )  , 

0  =  -f  (x)  +  g  (x )  +  C-l,' 


0  <  x  <  L, 


(22) 


Ik 


where  cx  is  an  integration  constant.  Adding  and  subtracting  these 
equations  now  gives 


f(x)  =  5*  , 

g(x)  , 


0  <  x  <  L  . 

—  9 


(23) 


Ci  evidently  contributes  nothing  to  (19),  and  will  hence  be  placed 
equal  to  zero,  giving 


f(x)  =  0,  g (x )  =  0  if  0  <  x  <  L  .  ( 2k ) 


Since  f(x-ct)  gives  a  wave  traveling  to  the  right,  and  g(x+ct) 
gives  one  traveling  to  the  left,  we  see  that  in  addition  to  the 
interval  0  <  x  <  L,  f(x)  must  be  determined  for  all  negative  values 
of  x,  and  g(x)  must  be  determined  for  all  positive  values  of  x. 
These  two  tasks  are  essentially  the  same,  however,  for  the  third 
condition  in  (20)  gives 


0  =  f(-ct)  +  g(ct)j 

(25) 

N'S) 

1 

II 

1 

Vi 

(26) 

The  last  condition  in  (20)  may  be  written 


s  =  vQt  when  x  =  L  .  (27) 

Equation  (19)  in  this  gives 

vQt  =  f(L-ct)  +  g(L+ct)  .  (28) 

Using  (26)  and  (28)  we  can  now  determine  f(x)  and  g(x)  as  follows. 


Let  us  construct  f(x)  and  g(x),  working  outward  on  both  sides 
from  the  interval  0  <  x  <  Lj  and  let  both  f(x)  and  g(x)  be  drawn 
using  the  same  axes.  Fig.  7.  We  now  proceed  as  follows. 
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Fig.  7  Functions  involved  in  the  general  solution  of  the  partial 
differential  equation. 

1.  Since  g(x)  =0  if  0  <  x  <  L,  f(x)  =  0  for  -L  <  x  <  0, 
due  to  (26). 

2.  Noting  from  (19)  that  s  is  given  in  the  interval  0  <  x  <  L 
by  letting  the  curve  f(x)  move  to  the  right  with  velocity  c,  and  the 
curve-  g(x)  move  to  the  left  with  velocity  c,  both  starting  at  t  =0; 
we  see  that  f  does  not  affect  s  at  x  =  L  until  it  has  moved  a  distance 
2L  to  the  right.  During  this  time  g  moves  a  distance  2L  to  the  left, 
and  in  order  to  maintain  (27)  must  slope  upward  with  a  slope  of  vQ/c 
while  passing  the  point  x  =  L.  g  must  therefore  slope  upward  with 
this  slope  in  the  interval  L  <  x  <  JL,  as  shown  in  Fig.  7- 

3.  It  follows  now  from  (26)  that  f  slopes  upward  with  slope 
vQ/c  in  the  interval  -3L  <  x  <  -L. 

4.  After  f  and  g  have  traveled  a  distance  2L,  f  arrives  at 

x  =  L  and  would  cancel  any  further  growth  of  s  at  this  point  if  the 
slope  of  g  were  not  altered.  In  order  to  maintain  the  rate  of  growth 
of  s  we  increase  the  slope  of  g  by  an  amount  vQ/c  to  cancel  the 
effect  of  f.  This  maintains  the  growth  of  s  at  x  =  L  while  the 
curves  move  an  additional  distance  2L;  g(x)  therefore  slopes  upward 
with  a  slope  2v0/c  in  the  interval  3L  <  x  <  5 L. 
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5.  It  now  follows  from  (26)  that  f  slopes  upward  with  a  slope 
2vq/c  in  the  interval  -  %  <  x  <  -  3L. 

6.  After  the  curves  have  moved  an  additional  distance  2L  the 
portion  of  f  with  slope  2v0/c  arrives  at  x  =  L,  and  it  becomes 
necessary  to  increase  the  slope  of  g  by  an  additional  amount  vQ/c 
to  give  a  slope  of  3vo/c,  and  so  on. 

We  now  see  that  g(x)  consists  of  a  broken  line.  Beyond  the 
first  interval  0  <  x  <  L,  in  which  g(x)  =0,  g(x)  consists  of  a 
series  of  straight  line  segments  of  horizontal  length  2L  and 
slopes  vQ/c,  2v0/c,  3vo/c>  kv0/c,  ...  respectively.  We  also  see  that 
f(x)  may  be  obtained  by  rotating  the  curve  of  g(x)  180°  about  the 
origin,  and  hence  consists  also  of  a  broken  line.  Fig.  7- 

To  obtain  s(x,t)  we  add  the  ordinate  of  g  at  abscissa  (x+ct) 
to  that  of  f  at  abscissa  (x-ct),  as  indicated  by  (19)*  We  see, 
however,  that  if  ct  >x,  the  ordinate  of  f  at  abscissa  x-ct  is  the 
negative  of  the  ordinate  of  g  at  abscissa  (ct-x).  Since  if  ct  <  x 
the  contribution  of  f  to  s  is  zero,  it  now  follows  that  f  (x)  can 
be  discarded,  and  s(x, t)  determined  from  g(x)  as  indicated  in 
Fig.  8. 
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If  we  reduce  ct  we  see  that  when  it  reaches  x  we  have  already 
lost  any  contribution  by  f.  The  above  construction  for  s  can  hence 
be  rendered  valid  for  values  of  ct  right  down  to  zero  by  extending 
g(x)  into  the  interval  -  L  <  x  <  0  and  placing 

g(x)=0  if  -  L  <  x  <  0.  (29) 


It  may  finally  be  mentioned  that  a  similar  procedure  could  be 
devised  if  vG  were  not  constant,  but  a  given  function  of  t. 

Fourier  Series  Solution. 

The  problem  we  have  just  solved  can  also  be  solved  by  using 
Fourier  technique.  The  result  obtained  will  be  less  general,  but 
it  will  be  needed  to  validate  some  of  the  practical  formulas  to  be 
obtained  later.  In  order  to  apply  Fourier  technique  it  is  necessary 
to  get  rid  of  the  nonhomogeneous  boundary  condition  at  the  right 
hand  end  of  the  string,  namely  the  condition 


^2  =  vn  when  x 

at  0 


=  L  . 


In  order  to  do  this  we  place 


v0  x  t 

s  =  -  +  s. 


and  consider  the  problem  of  determining  s1 . 
(17)  we  obtain 


62si  m  i_  agsi 

ax2  c2  a  t2  * 


(30) 

(31) 

Substituting  (31)  in 


(32) 


hence  sx  satisfies  the  same  partial  differential  equation  as  s. 
Writing  (31)  in  the  form 


s 


1  “ 


s- 


VoXt 

L 


(33) 


and  substituting  from  (20),  we  obtain  the  following  initial  and 
boundary  conditions  on  sx , 
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Sx  =  0  when  t  =  0 ,  0  <  x  <  L, 
ds.  vQx 

jr  = '  t  when  t  =  °> 

sx  =  0  when  x  =  0,  t  >  0, 

Si  =  0  when  x  =  L,  t  >  0. 


Proceeding  in  the  usual  manner  we  obtain  particular  solutions 
of  (32)  of  the  form 


sip  =  X(x)  T(t)  ;  (35) 

(35)  in  (32)  gives 

X"  T  =  i*  XT"  , 
c2 

or 

y»»  ^  m*1 

— =^2f-=-k2=  constant  (36) 

since  is  a  function  of  x  alone,  and  2—  is  a  function  of  t  alone. 

A  1 

Then  (36)  gives 


X"  +  k2  X  =  0  , 
T"  +  k2c2  T=  0 


} 


(37) 


Since  both  of  these  equations  are  linear  and  have  constant  coefficients 
we  can  write  down  the  solutions  immediately,  thus 

X  =  cx  sin  kx  +  cos  kx  , 

T  =  c3  sin  kct  +  c4  cos  kct  . 

Placing 

c2  =  0  >  1  .  . 

nw  7  (59) 

k  =  £2-  ,  n  =  1,  2,  3,  .  .  .  .  ,  J 

we  satisfy  the  end  conditions  in  (3k);  and  placing 

c4  =  0  (k0) 
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we  satisfy  the  initial  condition  on  displacement.  Superimposing  the 
particular  solutions  given  by  the  various  values  of  n  we  obtain 


v  „  nnx  nnct 

=  £  cn  sin  — = —  sin  ■  — 

n=l  L  L 


(hi) 


which  expression  satisfies  the  partial  differential  equation,  the  end 
conditions,  and  the  initial  condition  on  displacement.  To  satisfy 
the  initial  condition  on  velocity  we  place 


l  vat }. 


M 


or 


t=0 
nnc 


.  La*  -  r  c  £*£  sin 

L  'n=iCn  IT  Sln  T 

vQx  oo  .  nnx 

-  -2-  =  I  cnn  sin  — — 
nc  n=l  n  L 


nnx 


(*3) 


Hence  cnn  is  the  coefficient  of  sin  — in  the  Fourier  sine  series 

L 

expansion  of  (-vQx/nc).  It  follows  that 


2  f  vox  .  nnx  , 

-nn  =  t~  ~  -  sin  dx  , 

n  L  I  nc  L  ’ 

Jo 

r 

3  -  -  sin  25*  .  i*  cos  n£] 

n  Lnnc  \[njiJ  L  nn  LJ 


m 


or 


~n 


2v0L 

cos  nn 

n2n2c 

=  (-l)n  . 

n2n2c 


(^5) 


(*6) 


Placing  (U6)  in  (4l),  and  then  (4l)  in  (3l)>  we  obtain  finally 


s  =  v« 


(ii  +  SL  5  kil!  sin  sin  SSSi') 
V  L  fl2c  n=l  n2  L  L  / 


0*7) 


This  expression  gives  s  as  the  sum  of  a  steady  proportional  stretching 
and  a  set  of  vibrations  in  natural  modes  at  successively  higher  frequencies. 
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5.  Integral  Equation  Approach.  Conversion  of  the  Differential 

Equation,  and  Solution  by  Successive  Substitutions. 

Since  we  do  not  wish  to  limit  ourselves  to  the  solution  of  linear 
problems,  and  since  none  of  the  solutions  so  far  presented  furnishes 
a  suitable  means  of  dealing  with  nonlinear  or  time-dependent  materials, 
our  main  problem  is  not  yet  solved.  The  integral  equation  approach, 
which  will  now  be  presented,  has  in  the  past  led  to  the  solution  of  a 
variety  of  nonlinear  problems.  The  differential  equation  is  converted 
to  an  integral  equation,  and  the  integral  equation  is  then  solved  by 
successive  substitutions.  The  problems  solved  in  the  past  by  this 
technique  involved  but  a  single  independent  variable,  and  the  corres¬ 
ponding  differential  equations  were  ordinary;  whereas  in  the  present 
problem  two  independent  variables  are  involved,  and  the  differential 
equation  is  partial.  Nevertheless,  let  us  try  this  approach  and  see 
what  happens.  In  so  doing  it  is  desirable  first  to  devise  the  method 
and  apply  it  to  a  few  problems  whose  solutions  are  known,  before 
attacking  the  main  problem. 

Accordingly,  let  us  first  consider  the  case  where  the  string 
extends  from  minus  infinity  to  plus  infinity.  Noting  (17),  (l8), 
and  (12)  we  see  that  the  differential  equation  for  longitudinal 
motion  is 


(48) 


dt2  dx2 


where  for  the  present  we  consider  c  to  be  constant.  As  initial  con¬ 
ditions  let  us  place 


s(x,o)  =  f(x),  (rf)  =  0, 

8  “t=o 


(49) 


corresponding  to  which  the  initial  displacement  at  x  is  f(x),  and  the 
initial  velocity  is  zero.  Integrating  (48)  twice  with  respect  to  t 
we  obtain 


df-2  dt, 
dx2 


(50) 


s  =  s(x,o)  +  t(-^-)  +  c 

dtt=o 


2-|  dt  dt  . 
dx2 


(51) 
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With  the  initial  conditions  (49)  this  becomes 


s 


f(x) 


(52) 


Putting  s  equal  to  an  assumed  function  sQ  in  the  integrand  we  obtain 
a  function  Si  for  the  left  hand  side;  then  inserting  Si  in  the  inte¬ 
grand  we  obtain  s2,  and  so  on.  We  thus  have  a  method  of  successive 
substitutions  defined  by  the  following  set  of  equations 


Vi  ■ f(x) 


+  c 


if 


if 

dx2 


dt  dt,  n^^  1^  2^ 


(53) 


There  is,  of  course,  a  question  of  convergence.  We  hope  that 
ultimately  sn+-^  will  differ  inappreciably  from  sn  • 

In  order  to  start  the  process  let  us  place  sQ  =  0;  then  (53)  gives 


Sl 

=  f(x). 

s2 

=  f(x)  + 

c2t2 

T’- 

f"(x). 

S3 

=  f (x )  + 

catS 

2 

f"(x)  + 

C*t4 

TT 

S4 

=  f(x)  + 

• 

C2t2 

2 

f”(x)  + 

c*t4 

b: 

• 

sn 

=  f(x)  + 

eft2 

2.' 

f "(x)  + 

c4t4 

4: 

Ultimately 


f1V(x)  , 
f1V(x)  ♦ 


fiV(x)  + 


sn-s I  J 


s  •  f(x)  ♦  s|$2f(x)  +  slii 


cit4 

b: 


iv.  . 
f  (x) 


c6t6 


fVi(x) 


(55) 


But  one  form  of  Taylor's  series  is 

0  (x+h)  =  0  (x)  +  0'  (x)  h  +  0"  (x)fc|  +  0"'  (x)  £2  +  ...  .  (56) 

2.  3* 


I 


* 


* 
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It  follows  that 


f(x+ct)  =  f(x)  +  ct  f'(x)  +  f"(x)  +  ...  , 

<—  • 

f(x-ct)  =  f(x)  -  ct  f * (x )  +  s£if  f " (x )  -  ...  ; 

2 1 


} 


(57) 


hence  (55)  becomes 


-  IF 


2L 


f(x+ct)  +  f(x-ct) 


>]. 


(58) 


which  is  known  to  be  the  exact  solution  of  the  problem,  subject  to  the 
usual  conditions  of  convergence.  The  rate  of  convergence  of  the  series 
is  important;  if  the  rate  is  too  low  the  method  may  not  be  good  for 
practical  use. 

As  our  next  problem  let  us  consider  the  case  where  the  initial 
displacements  are  zero  and  the  initial  velocity  distribution  is  given, 
thus  instead  of  (49)  we  have 


s  (x,0)  -  0,  (&)  =  g(x). 


(59) 


t=o 


The  differential  equation  (48)  is  still  valid;  hence  integrating  twice 
with  respect  to  t,  as  before,  we  obtain  (51),  which  due  to  (59)  becomes 


s  =  t  g(x)  +  c' 


if. 


d- -  dt  dt . 
8x2 


(60) 


Proceeding  as  with  (55)  we  introduce  a  method  of  successive  substit¬ 
utions  described  by  the  equations 


5n+l 


*  g(x>  ♦  c’jT 


dx2 


dt  dt,  n  =  0,  1,  2,  ... 


(61) 


Starting  the  process  by  placing  sQ  =  0,  (6l)  gives 
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si  =  t  g(x). 


s2  =  t  g(x)  +  g"(x), 

J  • 


s3  =  t  g(x)  +  g"(x)  +  £lilglv(x) 

5*  5« 


*«+.5  iv, 


(62) 


,2*3 


sn  =t  g(x)  +  ^Ju  g"(x)  +  ...  + 

3  • 


2n-2,2n-i  2n-2,  N 
c_  t  (x) 


2n- 1. 


Ultimately 

S  =  t  g(x)  +  g"(x)  +  giv(x)  +  . . . 

Noting  (5 6)  and  (57)  we  see  that 


i; 

i: 


,2+3 


g  (x+ct)  dt  =  tg(x)  +  g'(x)  +  £_i_  g"(x)  +  ... 

2 :  5:  ' 


g(x-ct)  dt  =  tg (x )  -  ~  g'(x)  +  g"(x)  +  ...  . 

Ol  3!  > 
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hence  (63)  becomes 
s  -±/** 


Placing 


we  obtain 


r 


g(x+ct)  dt  = 


also,  placing 


(65) 


(64) 


(x+ct)  +  g(x-ct)J  dt  . 

(65) 

x  +  ct. 

> 

cdt. 

] 

^rtc+ct 

=  -J  g(u)  duj 

(67) 

A 

x-ct. 

) 

-cdt. 

|  (68) 
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we  obtain 


g(x-ct)  dt  =  - 


du. 


Putting  equations  (67)  and  (69)  in  (65)  gives  finally 


s 


x+ct 

g(u)  du 

-ct 


> 


which  is  known  to  he  the  exact  solution  of  the  problem. 


(69) 


(70) 


As  our  last  problem  of  this  kind,  let  us  consider  the  case  where 
the  initial  distributions  of  both  the  displacement  and  the  velocity 
are  given,  thus 


s(x,o)  =  f  (x). 


(— ) 


t»o 


=  g(x). 


(71) 


Equation  (51)  is  now  valid,  and  becomes 
s  =  f(x)  +  tg(x)  +  c2 


11 


dt  dt. 


(72) 


Our  method  of  successive  substitutions  is  now  described  by  the  equations 


sn+1  =  f(x)  +  tg(x)  +  c 


it* 


dt  dt,  n=0,  1,  2,  ...  .  (75) 


Starting  the  process  by  placing  sQ  =  0  we  obtain  for  sx,  s2, . . . , 
respectively,  the  sums  of  the  corresponding  expressions  given  in 
(54)  and  (62).  The  expression  ultimately  obtained  for  s  is  hence 
the  sum  of  those  given  in  (55)  and  (65);  and  this  can  evidently  be 
written  as  the  sum  of  (58)  and  (70).  We  thus  obtain  finally 


=  ijjf(x+ct)  +  f(x-ct)J  +  g (u )du 


(7M 


«/x-ct 


which  is  known  to  be  the  exact  solution  of  the  problem  for  the  linear 
case  [8]. 


25 


Extension  to  the  Case  of  a  Finite  Length. 


In  the  above,  in  order  to  fix  ideas,  the  piece  was  considered  to 
be  infinite  in  length;  however,  since  in  the  successive  substitution 
process  both  integrations  were  on  t  with  x  held  constant,  it  was  not 
necessary  to  consider  the  piece  to  be  infinite,  or  to  extend  outside 
of  the  interval  0  <  x  <  L.  Accordingly,  let  us  consider  the  problem 
as  it  is  actually  given.  That  is  the  piece  is  finite  in  length. 


extending  from  0  to  L;  the 

end 

conditions  are 

s(0, t )  =  0, 

<1 

)  =  v0  i 

x  =  L 

(75) 

and  the  initial  conditions 

are 

s(x,0)  =  0, 

0  <  x  <  L, 

(76) 

#)  -  0, 

dt  t=0 

0  <  x  <  L, 

(77) 

(bs\ 

fc’to  ’ T°! 

II 

X 

(78) 

bs 

We  note  that  (— 0  is  discontinuous  when  x  =  L  . 

dt  t=0 

The  integral  equation  is  given  by  (60)  or  by  (72)  with  f(x) 
placed  equal  to  zero  in  accordance  with  (76). 

We  note  that  g(x)  =  (—)  is  given  by  (77)  and  (78). 

ot  t=0 


Placing 


g(x)  =  +gi(x). 


(79) 


where 


gi(x)  if  0  <  x  <  L, 


gi (x)  =0  if  x  =  L, 


r 


(80) 


we  see  that  gi(x)  can  be  expressed  completely  —  including  the  discon¬ 
tinuity  at  x  =  L  —  by  the  sine  series 


gi(x)  =  Z  sin  — 
i— 1  L 


(81) 
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where 


*i  =  ei(x 


)  sin  dx. 


Substituting  from  (80)  this  becomes 

»L 


Jo 

vo  r rL  \2  .  inx  Lx  ijtx'“|L 

fa”  [<15  >  sin  T  -  15  cos  T Jo 


2v, 


i  2vr 


hence 


ai  =  Tn  cos  i5t  =  (-1)  1“  ; 


gi(x)  =  lli)  sin  . 

n  i=l  i  L 


This  in  (79)  gives 

00 

g(x)  kl  1  sin  1M  - 


*  i  =  1 


(82) 


(83) 

(84) 


(85) 


We  now  substitute  (85)  in  (60),  and  then  start  the  successive  sub¬ 
stitution  process  by  first  placing  s  =0.  It  may  be  noted  that  the 
first  term  in  (85)  at  no  stage  contributes  anything  to  the  integral- 
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which  agrees  with  (47 ) .  We  thus  see  that  the  nonhomogeneous  end  con¬ 
dition  at  x  =  L  can  be  used  without  alteration,  and  the  piece  con¬ 
sidered  to  extend  over  the  finite  interval  0  <  x  <  L. 

In  appendixes  A  and  B  a  method  is  provided  for  handling  problems 
in  which  the  motion  of  the  impacted  end  of  the  string  is  arbitrary. 
This  permits  us,  for  example,  to  specify  a  finite  velocity  rise-time, 
which  is  more  realistic  than  the  zero-rise-time  we  have  been  using. 
The  method  is  based  on  the  superposition  principle  and  requires  the 
material  to  be  linear.  The  development  is  given  in  Appendix  A  and 
the  corresponding  formulas  for  numerical  approximation  are  given  in 
Appendix  B. 

Since  the  method  of  successive  substitutions  has  been  shown  to 
be  valid  for  a  linear  (Hooke's  law)  string  of  finite  length,  we  can 
feel  confident  that  it  will  also  work  when  applied  to  nonlinear, 
time -dependent  materials.  The  problem  will  be  treated  directly  in 
the  form  in  which  it  is  given — that  is  with  a  piece  of  finite  length 
and  a  nonhomogeneous  boundary  condition.  It  is  likely  that  the  non¬ 
linearity  and  creep  will  not  greatly  affect  the  convergence  of  the 
process,  but  will  merely  complicate  the  numerical  work. 

6.  Development  of  Approximate  Methods  for  Solving  the  Integral 

Equation  by  Successive  Substitution. 

We  wish  to  obtain  the  solution  of  the  nonlinear  problem  with 
creep.  To  apply  the  theory  developed  above,  this  requires  the 
devising  of  approximate  methods  for  carrying  out  the  successive 
substitution  process.  In  so  doing  the  difficulties  encountered  will 
very  likely  be  about  the  same  regardless  of  whether  or  not  we  con¬ 
sider  nonlinearity  and  creepj  hence  in  order  to  keep  the  amount  of 
numerical  work  as  small  as  possible,  let  us  first  work  with  the 
linear  problem.  Since  we  have  the  solution  of  this  problem  we  can 
use  it  as  a  check  on  our  results,  which  is  an  additional  advantage. 

Nature  of  the  solution  of  the  Linear  Problem.  Referring  to 
Fig.  8  we  see  that  s(x, t)  is  given  by  the  increment  of  g  over  an 
interval  2x  whose  midpoint  is  at  abscissa  ct.  Since  g(x)  is  a 
broken  line  it  follows  immediately  that  in  general  ds/dx  and  ds/dt 
are  constant.  The  same  is  hence  true  of  the  vector 


i  + 
5x 


j 


5s 

dt 


-k, 


whose  direction  is  normal  to  the  surface  s(x,  t).  We  thus  see  that 
the  surface  s(x,t)  is  composed  of  pieces  of  planes. 


The  values  of  ds/dt  and  ds/dx  can,  in  addition,  be  determined. 
In  so  doing  we  note  that  there  are  two  possibilities,  (a)  and  (b). 


In  (a),  g(x)  has  no  corner  in  the  interval  of  width  2x.  In  this 
case  g(x)  has  a  constant  slope  nv0/c  in  the  interval.  Holding  x 
fixed  and  increasing  t  by  an  amount  dt  we  obtain 


ds  =  0  ; 

hence 

£§  =  0  . 
at 


(88) 


Next,  let  us  hold  t  fixed  and  increase  x  by  an  amount  dx;  then 

ds  =  2n  —  dx  ; 
c 


hence 


(89) 


Here  n  is  an  integer  which  is  determined  by  the  particular  branch  of 
g(x)  which  we  are  considering.  In  case  (b),  g(x)  has  a  corner  in  the 
interval  of  width  2x.  In  this  case  g(x)  consists  of  two  straight 
lines  in  the  interval,  their  slopes  being  n  v0/c  and  (n+l)  vQ/c, 
respectively.  Proceeding  as  before,  we  hold  x  constant;  then 


hence 


ds  =  (n+l)  -2  cdt  -  n  cdt 
c  c 

ds  =  vQ  dt; 


Next  we  hold  t  constant;  then 

ds  = 

ds  = 


(n+l)  —  dx  +  n  —  dx 
c  c 

(2n+l)  ^2  dx  ; 


(90) 
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hence 


|  -  <S»i>  •  (91) 

It  is  interesting  to  note  that  ds/dt  takes  only  the  values  0  and 
vQ,  whereas  ds/dx  keeps  increasing,  though  by  abrupt  increments  of 
vQ/c.  The  outer  boundary  of  the  surface  s(x,t)  consists  of  three 
straight  lines,  namely,  the  x  axis,  the  t  axis,  and  the  line 

s  =  v0t  ,  x  =  L.  (92) 

In  view  of  this  and  the  fact  that  planes  intersect  in  straight  lines 
we  see  that  the  various  plane  portions  of  the  surface  s(x,t)  are 
bounded  by  straight  lines. 

Referring  to  Fig.  8  we  see  that  as  t  increases  with  x  fixed  the 
transition  from  one  plane  portion  of  the  surface  to  the  next  occurs 
when  either  end  of  the  interval  of  width  2x  passes  a  corner  of  the 
broken  line  g(x).  Letting  x  differ  but  slightly  from  zero,  we  see 
that  the  transition  points  on  the  t  axis  occur  at 


or 


ct  — •  L,  JL,  % ,  ...  i 
*  _  L  2L  5L 

b  —  y  j  y  •  •  •  • 


(93) 


Similarly  if  we  let  x  be  slightly  smaller  than  L  we  see  that  the 
transition  points  on  the  boundary  line  (92)  are  at 


or 


ct  =  0,  2L,  kL,  ...  , 


t 


=  0, 


(9*0 


Applying  the  above  results,  noting  (93)  an<l  (9*0,  we  see  that  the 
transition  lines  between  the  various  plane  portions  of  the  surface 
s(x,  t)  are  as  shown  in  Fig.  These  lines  together  with  the  x  axis, 
the  t  axis,  and  the  line  (92)  form  the  boundaries  of  these  plane 
portions,  and  hence  determine  them.  A  three-dimensional  view  of  the 
surface  is  shown  in  Fig.  10. 


♦ 


•• 
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Fig.  9  Views  of  the  surface  s=s(x,t),  showing  transition  lines  between  the 
plane  portions  of  the  surface.  At  the  left  is  a  side  view  in  the 
+x  direction,  and  at  the  right  is  a  top  view  in  the  (-s)  direction. 
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Behavior  of  the  Integral  Equation.  On  each  of  the  plane  surfaces 
which  together  compose  the  surface  s(x,t)  we  have 


d2s 

aP  = 


0. 


(95) 


It  follows  that  in  the  integral  equation  (60),  namely 


s  =  tg(x)  +  c' 


u 


afs 

Bx2 


dt  dt. 


the  integrand  in  the  integral  vanishes  over  all  of  the  component  plane 
surfaces.  The  changes  in  the  first  of  the  two  integrals  hence  take 
place  abruptly  at  the  transition  lines  between  adjacent  surfaces. 

In  order  to  see  what  happens  at  such  a  line  let  us  consider  the 
cylindrical  surface 

z  =  f(x).  (9 6) 

Looking  down  from  above  let  us  consider  a  £  axis,  which  lies  in  the 
xy  plane  and  makes  an  angle  a  with  the  x  axis,  as  shown  in  Fig.  11. 


Fig.  11  Diagram  for  deriving  relations  satisfied  by  a 
cylindrical  surface. 
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Evidently 


x  =  ^  cos  a  ; 

(97) 

hence 

z  =  f(  J  cos  a). 

(98) 

It  follows  that 

=  f '  (  ?  cos  a)  cos  a, 

oj 

(99) 

=  f"(J  cos  a)  cos2  a  . 

(100) 

Similarly,  if  we  choose  an  T)  axis  perpendicular  to  the  J  axis. 

Fig.  11,  we  have 

x  =  T]  sin  a,  (101) 

=  f **"( T|  sin  a)  sin2  a  .  (102) 

dT 

Noting  (97)  and  (101)  we  see  from  (100)  and  (102)  that 

=  tl  tan2  a  .  (105) 

8n2  Zf 

Let  us  now  apply  this  result  to  the  first  transition  line,  which 
extends  from  point  x  =  L,  t  =  0  to  point  x  =  0,  t  =  L/c,  Fig.  12. 


Fig.  12  Application  of  cylindrical-surface  approximation  to  a  transition 
line. 


\ 


I 

« 
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Here  t  and  x  play  the  roles  of  f  and  q,  respectively;  and  we  consider 
the  slanting  plane  to  be  approximated  arbitrarily  closely  by  a 
cylindrical  surface  in  the  immediate  vicinity  of  the  transition  line, 
as  indicated  by  the  edge  view  Fig.  13. 


Fig.  13  Intersection  of  cylindrical 
-  surface  with  the  x,t  plane 


Equation  (IO3)  now  gives 

tan2  a.  (KA) 

dx2  dt2 

Considering  the  first  integration  in  the  double  integral  in 
(60),  we  have  seen  that  the  integral 


(105) 


receives  contributions  only  at  the  transition  lines.  It  is  hence  zero 
up  to  the  first  transition  line.  In  the  vicinity  of  this  line  we  sub¬ 


stitute  (104)  in  (105),  and  obtain 


tan2  a 


2 

g^dt.tan  C,A- 


(106) 


where  the  integral  extends  only  from  one  side  of  the  line  to  the  other. 
Here  A  (ds/dt)  is  the  increment  in  ds/dt  at  the  first  transition  line 
and  tan2  Ct  A(ds/dt)  is  the  corresponding  contribution  to  the  integral 
(105).  From  (90)  and  Fig.  10  we  see  that  this  increment  is 


also,  fl*om  Fig.  12  we  see  that 


:  vo  ; 

(107) 

=  - 

(108) 
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Substituting  in  (106)  we  obtain 


=  0  if  t  <  ti 
vo 

=  -5-  if  ti  <  t  <  t2 
c 


(109) 


where,  for  any  value  of  x,  tx  and  t2  are  the  values  of  t  for  the  first 
and  second  transition  lines,  respectively.  It  follows  that 


(110) 


Since  g(x)  =  0  if  x  <  L  it  follows  that  for  x  <  L  the  first  term  on 
the  right  hand  side  of  (60)  contributes  nothing,  s  being  given  entirely 
by  (110),  thus 


s  =  vQ  (t-ti). 


(Ill) 


Now  ti  depends  upon  x,  and  this  expression  is  evidently  in  accord  with 
the  results  shown  in  Fig.  10.  The  behavior  of  the  double  integral  at 
the  succeeding  transition  lines  can  be  obtained  in  a  similar  manner. 


Simplification  of  the  Integral  Equation.  Dividing  by  LvQ,  the 
integral  equation  (60)  may  be  written 


Placing 


(112)  becomes 


(112) 


(115) 


(111) 


But  this  integral  equation  is  identical  with  that  obtained  by  placing 


L  =  1,  c  =  1,  vQ  =  1 


(115) 


5  6 


in  (60).  It  follows  that  in  devising  approximate  methods  of  solution 
we  need  only  consider  this  special  case,  since  by  using  (113 )  the 
results  so  obtained  can  be  applied  to  any  case. 

Numerical  Integration  Scheme.  The  first  question  which  arises 
is  what  formulas  to  choose  for  numerical  differentiation  and  integ¬ 
ration.  In  order  to  keep  the  individual  steps  as  simple  as  possible 
let  us,  at  least  temporarily,  choose  the  trapezoidal  formula  for 
integration,  and  the  three  point  differentiation  formula  based  on 
the  quadratic  polynomial.  We  thus  have 

+  nh 

f(x)dx  =  h  [  -  f(a)  +  f(a+h)  +  f(a+2h)  +  ... 

2 


+  f(a+[n-l]h)  +  i  f(a+nh)], 

2 

(116) 

f”(x)  =h-2  [  f(x-h)  -  2f(x)  +  f (x+h)  ]  , 

(117) 

where  h  is  the  equal  spacing  between  adjacent  points.  Accuracy  will 
be  achieved  by  having  a  sufficiently  small  spacing  rather  than  by 
using  more  elaborate  formulas  for  differentiation  and  integration. 

Taking  advantage  of  (115)  we  shall  place  L  =  1,  c  =  1,  vQ  »  lj 

also 

g(x)  =  0  if  0  <  x  <  1, 
g(x)  =  1  if  x  =  1. 

Since  there  is  no  danger  of  ambiguity  we  shall  omit  the  bars  on  s, 
t,  x,  and  g;  however,  (113)  can  be  used  at  any  time  to  convert  any 
given  problem  to  that  which  we  are  now  considering.  Figure  9  now 
becomes  Fig.  14.  On  this  figure  the  values  of  ds/dx  and  ds/dt  on 
the  various  component  plane  surfaces  are  labeled,  these  values  being 
given  by  (88),  (89),  (90),  and  (91)*  For  convenience  we  shall  use 
the  same  spacing  h  in  both  x  and  t,  the  result  being  a  square 
lattice  in  the  xt  plane,  Fig.  14.  Later  some  preliminary  calcula¬ 
tions  will  be  made  with  h  =  .05;  however,  for  the  present  let  us 
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Fig.  l4  Top  view  of  the  s  (x,t)  surface  with 

variables  normalized.  Values  of  ds/dt 
and  ds/dx  are  shown  for  the  various  regions; 
the  values  change  abruptly  at  the  transition 
lines. 
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suppose  that  h  is  chosen  arbitrarily,  then  fixed.  On  the  boundaries 


s  =  0  on  the  x-axis  and  t-axis, 
s  =  t  on  the  line  x  =  1. 


} 


Since  g(x)  -  0  if  x  <  1  our  integral  equation  (114)  becomes 


s 


r4  dt  dt. 
bx 


(119) 


(120) 


In  the  integral  equation  (120)  we  see  from  (117)  that  the  inte¬ 
grand  contains  a  factor  1/h2.  however,  each  of  the  two  integrations 
using  (ll6)  introduces  a  factor  h,  the  result  being  that  h  disappears 
entirely  from  the  integral  equation.  We  are  therefore  justified  in 
replacing  h  by  1  in  our  calculations.  Of  course  n,  the  number  of 
subintervals  in  the  unit  length  of  string,  is  still  determined  by  the 
actual  value  of  h. 

In  carrying  out  each  repetition  of  the  successive  substitution 
process  it  would,  of  course,  be  possible  to  cover  the  entire  xt  plane 
0  <  x  <  1>  0  5  t  <  T>  where  T  is  the  greatest  value  of  t  in  which  we 
are  interested.  However,  since  for  the  larger  values  of  t  the  inter¬ 
mediate  solutions  may  fluctuate  wildly,  there  appears  to  be  no  point 
in  calculating  them  for  values  of  t  which  appreciably  exceed  that  for 
which  the  solution  has  converged.  In  fact  we  have  everything  to 
gain  and  nothing  to  lose  by  finishing  the  repetitions  for  one  value 
of  t  before  proceeding  to  the  next.  The  intermediate  solution 
obtained  is  hence  that  for  a  single  row  of  lattice  points  correspond¬ 
ing  to  a  constant  value  of  t.  Let  us  place 


sij  =  s(ih,  jh)  i  =  1,  2,  ...,  nj 

j  =  1,  2,  ...  $ 
n  =  1/h. 


(121) 


Then  the  successive  substitution  process  is  first  applied  to  deter¬ 
mine  sxi>  s£1,  •••>  1?  which  lie  in  the  same  horizontal  row; 

when  this  row  is  finished  it  is  applied  to  determine  s  ,  s  ,  ..., 
sn-i,2  >  an<*  so  on*  We  note  that  12  22 


=0  if  i  =  0  or  j  =  0, 


Snj  Jh* 


(122) 
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7 •  Preliminary  Investigations. 


\ 


Before  beginning  a  program  of  numerical  calculation  several  things 
were  done  by  way  of  preliminary  investigation  in  order  to  test  the 
convergence,  accuracy,  and  feasibility  of  the  successive  substitution 
process.  These  will  now  be  described  briefly.  Only  a  few  of  the  more 
significant  results  of  the  feasibility  studies  will  be  presented  in 
the  present  report.  All  of  the  results,  comprising  mathematical 
developments,  tables,  and  curves,  are  available  in  reference  [9]. 

Use  of  the  Exact  Solution  as  the  First  Approximation.  First,  as 
a  preliminary  check  the  exact  solution  s(x,t)  was  chosen  to  start  the 
successive  substitution  process,  and  the  result  of  one  cycle  of  this 
process  was  calculated.  Referring  to  the  oblique  lines  in  Fig.  14 
as  "transition  lines",  and  the  bands  composed  of  points  which  lie 
within  a  distance  hN2  of  these  lines  as  "transition  bands",  it  was 
shown  that  if  one  cycle  of  the  successive  substitution  process  be 
carried  out  starting  with  the  exact  solution,  the  approximate  solu¬ 

tion  obtained  differs  from  the  exact  solution  only  in  the  transition 
bands,  in  which  it  is  rounded  off  somewhat  when  compared  with  the 

sharp  corners  at  the  transition  lines  of  the  exact  solution.  The 

width  of  a  transition  band  is  hv2,  and  is  hence  a  measure  of  the 
spacing  in  the  lattice. 

It  was  thus  shown  that  one  cycle  of  the  successive  substitution 
process  does  not  lead  to  a  wild  result,  but  to  a  highly  satisfactory 
one.  Another  repetition  of  this  cycle  would  probably  result  in  a  » 

further  rounding  off  of  the  corners;  however,  this  investigation  was 
not  pursued  further. 

Determination  of  the  Solution  by  Difference  Equations.  In  the 
case  of  the  present  linear  problem  it  is  possible  to  determine  the 
solution  toward  which  the  successive  substitution  process  converges 
using  difference  equations.  The  solution  of  the  difference  equations 
gives  exactly  the  values  toward  which  the  successive  substitution 
process  converges,  provided  the  difference  equations  and  the  sub¬ 
stitution  process  are  both  based  on  the  same  lattice  spacing.  These 
common  values  differ  somewhat  from  the  exact  solution,  of  course. 

A  comparison  between  the  solution  of  the  difference  equations  and  the 

results  of  the  successive  substitution  process  is  practical  for  the  * 

first  few  rows  of  the  lattice.  The  comparison  is  made  one  row  at  a 

time. 


kO 


In  the  case  of  the  first  row  it  follows  from  (116),  (117),  (120), 
and  (122)  that  the  final  solution  satisfies  the  relations 


3 


*1-1,1  •  2  s 


’1-1,1 


6  s 


ii 


ii  +  Si+l,l)=  six 

(125) 

s .  =0. 

1+1,1 

(124) 

This  is  a  second  order  linear  difference  equation  with  constant 
coefficients,  which  was  solved  by  the  usual  procedure  for  handling 
such  equations. 

Not  only  the  quantities  s^x,  which  compose  the  first  row,  but  the 

quantities  s^j,  which  compose  the  j^*1  row,  can  be  determined  by  a 
difference  equation.  This  equation  is 

8i+i, J  "  6  8ij  + 

J-l 

=  “  ^  (8i+a  k  ‘  2  sik  +  si-x  k)>  (125) 

k=l 

which  is  also  a  second  order  linear  difference  equation  with  constant 
coefficients.  Quantities  pertaining  to  the  first  (j-l)  rows  are  at 
this  point  known.  This  equation  was  also  solved  by  the  usual  methods. 
We  note  that  the  above  difference  equations  have  nothing  in  common 
with  those  obtained  by  replacing  the  partial  derivatives  in  (48)  by 
their  difference  approximations. 


The  values  of  s^j  for  the  first  six  rows  were  calculated  by 
means  of  the  above  difference  equations.  The  existence  of  a  solution 
of  the  integral  equation  (120)  was  thereby  established;  furthermore, 
the  values  of  s^  so  obtained  were  available  as  a  standard  of  compari¬ 
son  to  evaluate  the  results  obtained  at  any  stage  of  the  successive 
approximation  procedure — at  least  for  j  =  1,2,  ...,  6.  In  addition, 
we  could  see  whether  or  not  the  result  of  the  successive  substitution 
process,  when  sufficiently  converged,  gave  a  reasonably  good  approxi¬ 
mation  to  the  exact  solution.  From  Fig.  15  we  see  that  it  does; 
however,  we  should  expect  that  a  much  better  approximation  would  be 
obtained  if  the  motion  of  the  moving  end  did  not  have  a  discontin¬ 
uity  in  velocity,  but  were  more  realistic.  As  things  are,  the 


X 

l 


Fig.  15  Comparison  between  the  approximate  (successive  substitution)  solution  and  the  exact 
solution,  for  n  =  20  and  j  =  1  to  6  (t  =  .05  to  .25).  The  approximate  solution 
agreed  exactly  with  the  solution  of  the  difference  equations,  which  showed  that 
the  successive  substitution  process  had  converged. 


approximate  solution  is  trying  to  follow  the  corners  which  appear 
in  the  graph  of  the  exact  solution. 

8.  Preliminary  Calculations. 

Various  preliminary  calculations  were  carried  out  as  follows,  using 
a  ten  bank  desk  calculator. 

A.  Calculation  of  by  Successive  Substitutions.  Sixteen  repeti¬ 
tions  of  the  successive  substitution  process  for  determining  B^1, 
the  values  of  Sj_ j  for  the  first  row,  were  carried  out.  The  results 

obtained  for  i  =  19  and  18  are  plotted  as  solid  lines  in  Fig.  16. 

Note  that  n  =  20  and  that  the  plotted  results  refer  to  the  two  lattice 
points  immediately  adjacent  to  the  moving  end  of  the  string.  The 
values  plotted  are  those  resulting  from  the  (k-l)th  repetition,  these 
being  used  to  begin  the  k1^  repetition.  The  dotted  and  dashed  lines 
in  the  figure  refer  to  modifications  of  the  successive  substitution 
method,  made  in  attempts  to  speed  up  convergence.  These  will  be  dis¬ 
cussed  later. 

In  determining  the  displacements  (si;L)  in  the  first  row,  we  may 
expect  convergence  first  for  the  higher  values  of  i  such  as  i  =  19 
or  18  that  are  near  the  given  value  S2Q  , .  Only  after  convergence 
in  this  region  is  well  on  the  road  can  we  expect  convergence  for  lower 
values  of  i;  and  such,  indeed,  is  the  case. 

The  values  of  s^1  obtained  after  sixteen  repetitions  were  com¬ 
pared  with  the  exact  values,  which  are  known.  It  appears  that  the 
process  is  converging;  however,  at  this  point  values  of  for 
i  <  17  have  not  yet  really  begun  to  settle  down. 

B.  Application  of  the  Successive  Substitution  Process  to  the  Results 

Obtained  for  the  First  Row  (.1=1)  by  the  Difference  Equation  Method. 

Starting  with  the  values  s^1  obtained  by  the  difference  equation 
method  for  i  =  9,  10,  ...,  19,  one  cycle  of  the  successive  substitu¬ 
tion  process  was  carried  out  to  obtain  values  of  si:L  for  i  =  10,  11, 

...,  19.  These  were  compared  with  the  original  values,  and  it  was 
found  that  to  the  capacity  of  the  machine  the  output  of  the  process 
duplicates  the  input.  This  provides  an  independent  check  on  the  values 
obtained  by  the  difference  equation  method;  also  it  illustrates  the 


Fig.  16  Convergence  of  s-values  at  two  selected  lattice  points  in  the  first  row.  Solid  line, 
unmodified  procedure}  other  lines,  modified  output  taken  as  a  x  input  +  0  unmodified 
output . 


last  step  which  would  have  been  taken  in  Part  A  if  the  successive 
substitution  process  had  been  carried  to  completion. 

C.  Calculation  of  s^  Using  a  Weighting  Factor  of  .3  with  the  Input 

and  .7  with  the  Output.  As  we  examine  Fig.  16  we  note  that  the  con¬ 
vergence  obtained  is  of  the  nature  of  a  damped  oscillation  -  the 
successive  calculated  values  overshoot  the  true  value.  In  view  of 
this  it  appears  that  the  true  value  would  be  more  closely  approxi¬ 
mated  by  an  average  of  the  values  obtained  after  two  consecutive 
repetitions  of  the  successive  substitution  process  than  by  the  value 
obtained  last.  The  procedure  used  so  far  may  be  shown  schematically 
by  the  diagram  Fig.  17* 
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Fig.  17  Schematic  representation  of  successive  substitution  process. 


Here  the  blocks  represent  the  various  repetitions^ 
the  kth  repetition  and  is  the  output  of  the  kth 
the  process  which  we  have  been  using 


1^  is  the  input  to 
repetition.  In 


To  increase  the  damping  of  the  oscillation  the  output  of  one 
stage  can  be  coupled  to  that  of  the  preceding  stage,  as  indicated 
by  Fig.  18,  and  described  by  the  relation 


a  J 

k-i 


+  S  J 

k 


=  I 

k+i 


(126) 


where  the  weighting  factors  a  and  P  are  constant,  .are  chosen  once  and 
for  all,  and  satisfy  the  relation 

a  +  p  =  1. 

^5 


(127) 
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Fig.  18  Schematic  representation  of  one  step  in  modified 
substitution  process,  in  which  the  weighting 
factors  a  and  p  are  used. 

In  choosing  a  and  0  a  number  of  factors  are  evident  as  follows. 

1.  If  a  =  1,  P  =  0  the  repetition  has  no  effect;  hence 
we  require  that  a  <  1,  p  >  0. 

2.  If  convergence  is  achieved  through  a  damped  oscillation 
when  a  =  0,  p  =  1,  the  damping  or  rate  of  convergence  will 
be  improved  if  P  is  reduced  from  unity  to  some  value  which 
exceeds  .5  (see  Fig.  1 6). 

3.  If  convergence  is  achieved  without  oscillation  by  approach¬ 
ing  the  limit  from  one  side  when  a  =  0,  P  =  1,  the  rate  of 
convergence  will  be  improved  if  P  is  increased  from  unity 
to  as  large  a  value  as  can  be  used  without  causing  over¬ 
shooting  and  hence  oscillation. 

k.  Since  the  requirements  of  2.  and  3*  are  contradictory,  and 
since  in  any  case  we  do  not  know  whether  or  not  the  con¬ 
vergence  is  achieved  through  a  damped  oscillation,  we  should 
choose  a  value  of  between  .5  and  1  which  is  small  enough  to 
greatly  help  the  convergence  if  oscillation  is  involved, 
but  not  so  small  as  to  greatly  hurt  the  convergence  if 
oscillation  is  not  involved.  Of  the  possible  choices 
p  =  .5,  .6,  .7,  .8,  .9,  we  can  rule  out  .5  as  being 
smaller  than  we  need  or  want,  and  .9  as  not  differing 
enough  from  unity.  If  there  is  oscillation,  P  =  .6 
is  very  likely  a  good  value,  but  .7  should  be  almost 
as  good  when  oscillation  is  present  and  appreciably  better 
if  there  is  no  oscillation.  We  could  put  P  =  .8 
but  this  seems  a  little  close  to  unity.  All  in  all,  it  ap¬ 
pears  that  the  best  value  we  can  choose  in  the  absence  of  any 


data  is  0  =  .7*  We  shall  hence  tentatively  place 


a  =  .3,  0  =  .7 


(128) 


We  now  have  two  ways  in  which  intermediate  results  are  improved, 
as  follows. 

a)  The  carrying  out  of  a  repetition  of  the  successive  sub¬ 
stitution  process. 

b)  The  calculation  of  the  weighted  average  (126)  to  give 
Ik+i* 

At  this  point  Robert  Crout*  suggested  that  since  the  weighted  average 
Ik  is  better  than  the  preceding  output  the  latter  should  be 

discarded  in  favor  of  the  former  in  subsequent  calculations,  in 
particular,  this  should  apply  to  (126).  This  suggestion  seems 
reasonable;  hence  we  shall  supersede  (126)  by  the  relation 


(129) 


The  considerations  which  led  to  the  choice  a  =  .3,  (3  =  .7  remain  un¬ 
altered;  hence  we  shall  tentatively  retain  (128). 

Since  the  calculation  of  the  weighted  average  b)  is  much  easier 
than  the  carrying  out  of  a  repetition  a),  and  does  lead  to  an  improve¬ 
ment  in  the  results,  it  was  suggested  by  Robert  Crout  that  where  the 
process  of  successive  repetitions  is  terminated  before  convergence 
has  been  achieved  the  last  step  be  of  the  type  b).  This  is  equivalent 
to  regarding  the  result  obtained  after  the  kt  repetition  to  be  1^ 

instead  of  J>.  This  suggestion  seemed  reasonable,  and  was  followed. 


Equations  (128)  and  (129)  together  give 
*5Ik  +  .7Jk  =  Ik+1- 


(130) 


Robert  Crout  made  many  of  the  preliminary  calculations  of  the 


present  investigation  on  a  desk  calculator,  under  the  direction 
of  his  father,  Prescott  D.  Crout. 


Fourteen  successive  repetitions  based  on  (I30)  were  carried  out. 

The  results  after  the  various  repetitions  are  shown  for  s18  x  and 
si8,i  in  Fig*  16.  We  see  that  the  present  method  of  successive 
repetitions  of  the  process  based  on  (130)  is  very  much  better  and 
leads  to  much  more  rapid  convergence  than  does  the  unmodified 
successive  substitution  process  considered  in  Part  A. 

Finally,  it  may  be  mentioned  that  we  should  expect  to  obtain 
convergence  using  any  value  of  3  lying  in  the  range  0  <  f3  <  1; 
however,  as  we  have  seen,  the  rate  of  convergence  depends  upon  the 
value  chosen. 

D.  Calculation  of  s„  Using  a  Weighting  Factor  of  .3  with  the  Input 
and  .3  with  the  Output.  Fourteen  repetitions  were  carried  out  in 
computation  of  s. ^  with  a  =  .5,  3  =  .5*  The  resulting  convergence 
of  s19>i  and  s18, 1  are  shown  in  Fig.  1 6.  We  note  that  the  convergence 
is  not  oscillatory.  We  have  seen  that  reducing  3  increases  the 
damping  of  an  oscillation. 


Comparing  the  results  of  parts  C  and  D  we  see  that  the  convergence 
is  better  with  a  =  .3,  3  =  *7  than  it  is  for  a  =  .5,  3  =  The  latter 
values  are  little,  if  any,  better  than  the  former  in  the  case  of 
oscillation,  and  are  appreciably  worse  if  there  is  no  oscillation. 


E. 


Calculation  of  s. 

-16 


Using  Successive  Substitutions. 


We  have  just 

finished  considering  the  determination  of  the  first  row  s^.  There 
is  no  reason  to  make  preliminary  calculations  on  more  than  one  addit¬ 
ional  row;  hence  let  us  proceed  directly  to  the  determination  of  the 
sixth  row,  since  that  is  the  last  for  which  we  have  calculated  exact 
values  of  s^.,  in  this  case  s^g,  using  the  difference  equation  method. 


We  shall  first  use  the  method  of  successive  substitutions  (a  = 
0,  3  =  l)>  end  shall  consider  the  exact  values  of  s^  for  the  first 
five  rows  to  be  given.  To  start  the  process,  values  of  s^5  in  the 
fifth  row,  which  are  known,  are  taken  as  the  first  approximation  to 
the  corresponding  values  of  siQ.  Fourteen  repetitions  were  carried 
out,  the  resulting  convergence  of  s19  6  and  s18  6  being  shown  in 
Fig.  19  (solid  lines).  In  this  figure  it  should  be  noted  that  the 
origin  is  not  shown.  The  repetition  number  is  (k-1). 
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Convergence  of  s-values  at  two  selected  lattice  points  in  the  sixth  row.  The  use  of  an 
extrapolated  input  is  compared  with  the  use  of  the  values  in  the  preceding  row.  Also, 
the  effect  of  using  weighting  factors  is  shown. 


Fig.  19 


F .  Calculation  of  g  Using  Successive  Substitutions  Starting  with 
Values  Obtained  by  Extrapolation.  The  present  calculations  differ 
from  those  in  Part  E  only  in  that  more  accurate  initial  values  are  used 
to  start  the  successive  substitution  process-  Whereas  in  Part  E  s.,5 
is  taken  as  the  first  approximation  to  s^  ,  we  shall  now  take  as  the 
first  approximation  to  s^6  the  value  obtained  by  linear  extrapolation 
using  s^4  and  s^g.  We  thus  have  as  our  first  approximation 

Si6  *  Si5  +  (Si5-Si4)  “  28i5-Si4- 

Four  repetitions  are  carried  out,  the  resulting  convergence  of 
S19  q  and  sie,  e  being  shown  in  Fig.  19*  The  repetition  number  is 
(k-1). 


Comparing  these  last  two  results,  it  appears  that  whatever  ad¬ 
vantage  is  obtained  by  using  the  more  accurate  initial  values  given 
by  (ljl)  can  be  achieved  more  easily  by  using  the  values  of  si5  as 
first  approximations  to  the  values  of  s.  and  using  two  or  three 
more  repetitions  of  the  iterative  process. 

G.  Calculation  of  s^  Using  a  Weighting  Factor  of  .3  with  the  Input 
and  -7  with  the  Output,  and  Starting  with  Values  Obtained  by 

Extrapolation.  Four  repetitions  were  carried  out  using  a  weight¬ 
ing  factor  of  .3  with  the  input  and  .7  with  the  output  (a  =  .3,  0  =  -7) • 
The  first  approximation  was  obtained  by  linear  extrapolation  using  (131). 
We  have  already  seen  in  section  F  that  the  linear  extrapolation  is  not 
worth  the  trouble  but  now  we  wish  to  see  the  advantage  of  using  various 
weighting  factors.  The  convergence  of  s19^6  and  s18^6  with  a  =  .3  and 
3  =  .7  is  shown  in  Fig.  19*  We  see  that  the  convergence  is  very  much 
improved  by  the  use  of  weighting  factors,  as  compared  with  the  con¬ 
vergence  obtained  in  Parts  E  and  F  using  the  method  of  successive  sub¬ 
stitutions  (a  =  0,  3  =  l). 

H.  Calculation  of  s^.  Using  a  Weighting  Factor  of  .3  with  the  Input 
and  .3  with  the  Output,  and  Starting  with  Values  Obtained  by 

Extrapolation.  Four  repetitions  were  carried  out  using  a  weighting 
factor  of  -5  with  the  input  and  .5  with  the  output  (a  =  .5,  0  =  »5)- 
The  first  approximation  was  obtained  as  in  F  and  G  by  linear  extrap¬ 
olation  using  (131).  The  resulting  convergence  of  s19^6  and  s18^6  is 
shown  in  Fig.  19 .  Comparing  the  results  in  Fig.  19  we  see  that  for 
S19, s  and  s18>6  there  is  not  much  difference  in  convergence  rates 


when  Ct  =  .3>  P  -  «7  and  when  Ct  =  .5,  0  =  .5-  However,  both  sets  of 
weighting  factors  give  much  more  rapid  convergence  than  the  unmodified 
process  that  corresponds  to  a  =  0,  0  =  1,  and  we  recall  that  s19^  x  and 
si8,i  converged  most  rapidly  when  a  =  .3,  3  =  .7.  As  a  consequence  of 
the  above  preliminary  calculations  we  have  found  that: 

1.  The  convergence  of  the  successive  substitution  process  is  for 
all  practical  purposes  established. 

2.  The  rate  of  convergence  is  very  greatly  improved  by  using 
weighting  factors,  as  described  on  page  U7. 

3.  On  the  basis  of  the  work  done  so  far  the  recommended  values 
of  the  weighting  factors  are  a  =  .3  and  0  =  .7,  which  lead 
to  the  fundamental  relation  (130). 

h.  The  use  of  linear  extrapolation  using  (131)  to  obtain  the 
first  approximation  is  not  worthwhile,  since  the  improve¬ 
ment  which  it  gives  can  more  easily  be  achieved  by  a  few 
more  repetitions.  It  is  recommended  that  the  known  values 
s^j^  be  taken  as  the  first  approximation  of  the  corres¬ 
ponding  values  Sjj. 

The  above  recommendations  are  incorporated  in  the  following  set 
of  instructions  for  programming  the  successive  repetition  process 
for  use  in  connection  with  a  digital  computer.  In  these  instructions 
A  is  a  measure  of  the  first  time- integral,  whereas  B  is  a  measure  of 
the  second. 

9.  Description  of  the  Repetitive  Process  When  Applied  to  the  Linear 

String. 

Our  problem  is  to  determine  the  values  of  the  quantity  s^  at  the 
points  of  a  square  lattice,  in  which  i  takes  the  values  0,  1,  2,  3,..., 
n;  and  j  takes  the  values  0,  1,  2,  ....  SLnce  n  is  given,  the  range 

of  values  taken  by  i  is  bounded.  The  range  of  values  taken  by  j, 
however,  has  no  upper  bound  except,  possibly,  one  which  will  be 
assigned  later  from  practical  considerations.  Fig.  20. 

The  values  of  sij  are  given  for  those  points  which  lie  on  the 
border,  thus 

sQj  =  0  j  =  0,  1,  2,  3,  ...; 

sio  =  0  1  =  °>  2,  3,  .  •  • ,  nj 

snj  is  given'  J  =  °>  2>  3,  •••  • 
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Fig.  20  Lattice  of  points  at  which  calculations  are  to  be 
made. 


The  values  of  s  _  at  the  other  points  are  determined  a  row  at  a  time  — 
that  is  for  j  fixed  and  1*1,  2,  3,  ...  ,  (n-l).  The  values  of  j 
are  taken  in  the  order  1,  2,  3,  ...  .  In  view  of  this  it  will 
evidently  be  sufficient  to  describe  how  to  obtain  the  values  of 
sij  for  the  row  when  the  values  for  all  of  the  preceding  rows 
are  known. 

"th 

The  values  of  s  for  the  j  row  are  obtained  as  the  limits 
gotten,  respectively,  by  applying  repeatedly  a  process  shortly  to 
be  described.  At  the  beginning  of  the  k^*1  repetition  we  have  the 
quantities 

(k)  (k)  (k)  (k) 

sij'  Ssi’  Bal’  ”*  '  sn-i, j  ^53) 


% 


i 


P 
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as  approximations  to  the  quantities 


V 


(13*0 


respectively.  Here  (k)  is  a  superscript,  not  an  exponent.  When  the 
process  is  applied  starting  with  the  quantities  (133 ),  the  quantities 


(k+i)  (k+i)  (k+i)  (k+i) 

SU  '  Saj  '  Saj  '  '*•  '  (135) 

are  obtained  as  approximations  to  (13*0  respectively.  The  sequence 
of  repetitions  is  started  by  taking  the  values  of  s^ .  for  the 
(j_l)  n  row  as  approximations  to  the  values  for  the  Jjt  row,  res¬ 
pectively,  thus 


lj  1*J-1  2  j 


y 


s 


(1) 

n-i, 


J 


Sn-i,J-x* 


(136) 


These  values  are,  of  course,  known.  The  sequence  of  repetitions  is 

terminated  when,  for  sufficiently  large  k,  the  output  (135)  duplicates 

the  input  (133)  bo  the  required  accuracy.  The  common  values  which 

then  compose  (133)  and  (135)  are  taken  as  the  final  result  (13*+)  •  The 

values  of  s.  .  which  compose  the  Jth  row  are  thus  determined. 

J 


Finally,  we  shall  describe  the  process  by  which  we  start  with 
(133)  and  obtain  (135) •  During  the  course  of  the  calculations  we 
compute  and  record,  in  addition  to  the  quantities  s„,  the  auxiliary 
quantities  and  B^.  These  quantities  do  not  exist  at  the  points 
of  the  left  border  of  the  lattice,  for  which  i  =0,  or  the  right 
border,  for  which  i  =  n.  Along  the  lower  border,  where  j  =  0,  they 
are  defined  to  be  zero,  thus 


A  =  0,  A  =  0,  A  =0,  . . .,  A 

10  ’  20  *  30  ’  ’  n-i 


=  0 


B  =  0,  B  =  0,  B  =  0,  . . .  ,  B 
10  20  30  u  1,0 


=  0. 


(137) 

(138) 


As  is  the  case  with  s.  the  auxiliary  quantities  A  and  B,  ,  are 

obtained  one  row  at  a  time,  are  known  in  the  first  (j-1;  rows,  and 
a^e ^obtaii^ecj.  in  the  j"fch  row  as  the  limits  approached  by  the  quantities 

A  .  and  B.,  .,  respectively,  during  the  course  of  the  successive 
i  j  J 

repetitions.  The  desired  basic  process  is  now  described  by  the 
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relations 


.00  .  00 

AiJ  '  Ai#J-i  ~si,J-i  "sij 


2  ^Si-i#j-i  +  Si+i,j-i 


+  s 


oo 

i-i#j  T  °i+i,j 


(k)  \ 

+  s,  .  4), 


i  -  1>  2#  • • ■ ;  (n-l)j 
i  =  1^  2,  •  •»#  (n-1); 


(139) 

(140) 


(k+i) 

8U 


•  7 


1=1#  2#  •  *•#  (n-1) . 


(141) 


Starting  with  the  given  quantities  s^  ,  j  being  fixed  and  i  taking  the 
values  1,  2,  ...  ,  (n-1)#  we  compute  the  quantities  Ajk'  using  (139)$ 
then  the  quantities  B^k  using  (140);  and  finally  the  quantities  sjk+1^ 
using  (l4l).  This  process  is  now  repeated  until  all  values  of  s^  have 
adequately  converged.  Row  j+1  is  then  computed  in  a  similar  manner  and 
the  calculation  of  successive  rows  is  continued  until  the  time  range  of 
the  problem  has  been  covered. 


10.  Typical  Results  Obtained  by  Digital  Computer,  Using  the  New  Methods. 

The  preceding  analysis  has  dealt  with  the  problem  of  a  linear  string# 
one  end  of  which  is  instantaneously  set  in  motion  with  constant  velocity 
at  time  0.  A  computer  program  applying  the  successive- substitution 
method  to  this  problem  was  written  and  tested  out.  Putting  n  =  20#  the 
computed  results  should  agree  exactly  with  the  solution  of  the  finite- 
difference  equations  shown  in  Fig.  15;  this#  in  fact#  they  did.  But 
whereas  with  the  finite  difference  equations  and  the  desk  calculator 
it  was  impractical  to  go  much  beyond  n  =  20  and  j  =6#  the  high-speed 
computer  could  easily  handle  n  =  100#  with  j  ranging  up  to  J>00  or  even 
higher . 

Figure  21  shows  displacement  s  versus  x  as  obtained  by  computer# 
with  n  =  100,  for  two  typical  rows  of  computer  output#  for  which 
j  =  75  and  J  =  275#  respectively.  The  computed  quantities  are 
normalized  according  to  equations  (114)  to  (117);  unit  strain  travels 
with  unit  velocity  in  a  string  of  unit  length.  Hence  the  first  curve 
(j  =  75)#  corresponds  to  t  =  .75#  when  the  wave  front  has  traveled 
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3 A  the  length  of  the  string.  The  second  curve  (j  =  275)  corresponds 
to  t  =  2.75>  when  the  wave  front  has  been  reflected  twice  and  has 
completed  3 >A  of  its  third  passage.  The  exact  solution  of  this 
problem  is  known  to  consist  of  two  straight  lines  having  a  corner 
(discontinuity  in  slope)  at  the  position  of  the  wave  front.  The 
differences  between  the  computer  solution  and  the  exact  solution  are 
hardly  noticeable  although  on  a  large-scale  graph  it  can  be  seen  that 
the  computer  solution  oscillates  from  side  to  side  of  the  exact 
solution. 

In  addition  to  the  displacement  s  it  is  almost  always  desirable 
to  compute  the  strain  e  =  ds/dx.  Differences  between  experimentally 
observed  and  calculated  strains  are  more  convenient  to  compare  than 
the  corresponding  differences  in  displacements.  Although  appearing 
in  the  basic  differential  and  integral  equations,  the  strain  does 
not  appear  explicitly  in  the  program  for  the  linear  string  where 
approximate  formulas  are  used.  Rather  the  quantities  from  which 
the  strain  could  be  calculated  are  to  be  found.  However,  in  later 
programs  various  quantities  are  functions  of  strain  and  also  strain 
values  are  required  for  comparison  with  other  calculations  and 
experimental  results.  Therefore  the  strain  results  for  this 
problem  should  be  examined.  They  were  calculated  by  hand  from  a 
simple  two  point  formula, 

h -i,3  *  "  (%J  '  >’ 

(Practically  the  same  results  were  obtained  with  a  five  point  formula. 

Figure  22  shows  strain  plotted  against  x  for  the  same  times  as 
in  Fig.  21.  The  exact  solution  is  shown  by  the  horizontal  and  verti¬ 
cal  straight  lines.  In  this  figure  the  departures  of  the  computer 
solution  from  the  exact  solution  are  very  much  more  noticeable  than 
they  were  in  Fig.  21.  This  is  to  be  expected,  since  the  ordinates 
in  Fig.  22  are  the  slopes  of  the  corresponding  curves  in  Fig.  21. 

The  oscillations  of  the  computer  solution  undoubtedly  originate 
because  a  finite -difference  method  cannot  conform  exactly  to  the 
sudden  discontinuity  in  strain  at  the  wave  front.  In  the  computer 
solution  the  wave  front  has  a  finite  slope  whereas  the  true  slope 
is  infinite.  In  the  computer  solution  the  strain  does  not  reach 
its  full  value  when  the  exact  solution  does,  but  lags  behind.  When 
the  computer  solution  does  reach  the  true  value  it  overshoots.  The 
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Fig.  22  Values  of  strain  (e )  taken  from  the  computer  data  of  Fig.  21. 
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overshooting  is  eventually  corrected  but  the  correction  is  too  large 
and  causes  undershooting.  The  result  is  that  the  computer  solution 
oscillates  back  and  forth  across  the  line  representing  the  exact 
solution. 

The  oscillations  of  the  computer  solution  above  and  below  the 
exact  solution  in  Fig.  22  are  larger  than  we  would  like.  Therefore, 
it  is  an  important  point  that  in  practical  problems  factors  are 
present  that  will  reduce  the  amplitude  of  the  oscillations.  To  see 
that  nonlinearity  of  the  stress- strain  curve  can  be  beneficial  in 
this  respect,  we  note  that  in  a  non-linear  material  the  velocity 
of  strain  propagation  depends  on  the  level  of  strain.  In  the  usual 
case  small  strains  travel  more  rapidly  than  large  strains.  In  such 
a  case  a  large  strain  consists  of  a  rapidly  moving  leading  edge 
followed  by  a  slower-moving  crest;  the  strain  rises  gradually  as 
the  wave -front  passes,  not  discontinuously.  The  slope  of  the  front 
becomes  more  and  more  gradual  as  the  wave  travels,  and  solutions 
computed  by  the  method  of  successive  substitutions  can  conform  more 
and  more  closely  to  it.  If  a  square  wave  were  ever  present  in  such 
a  case  it  would  rapidly  deform  into  a  wave  with  a  slanting  front. 

Situations  can  arise,  in  theory  at  least,  (when  the  stress- 
strain  curve  is  concave  upward)  in  which  large  strains  travel  more 
rapidly  than  small  strains.  In  such  a  case  an  initially  slanting 
wave  front  becomes  steeper  until  a  shock  wave  is  formed.  Such 
situations  will  be  encountered  later,  but  only  in  the  computations, 
not  in  the  experimental  data.  So  far  as  we  are  aware,  there  is  no 
conclusive  evidence  that  shock  waves  are  ever  produced  in  tensile 
impact  experiments.  Situations  in  which  the  computations  lead  to 
shock  waves  must  be  excluded  when  we  say  that  nonlinearity  reduces 
the  oscillations  in  computer  solutions;  some  situations  to  be 
treated  in  a  later  report  indicate  that  the  tendency  to  form  shock 
waves  may  aggravate  the  oscillations  in  computer  solutions. 

The  effect  of  time -dependence  (creep  or  recovery)  in  reducing 
oscillations  in  computer  solutions  should  always  be  favorable. 

When  part  of  the  response  of  the  material  is  delayed,  the  shape  of 
a  square  wave  will  be  distorted  as  the  wave  travels,  and  the  dis¬ 
tortion  will  tend  to  smooth  out  discontinuities,  making  it  easier 
for  the  successive- substitution  method  to  follow  the  actual  strains. 

In  the  exact  solutions  presented  above,  which  we  wish  to  match 
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with  the  computer  solution,  the  acceleration  of  the  moving  end  of  the 
string  occurs  in  zero  time.  If  we  give  up  this  condition,  which  is 
actually  unrealistic,  and  allow  the  end  of  the  string  to  experience 
a  constant  acceleration  during  a  small  finite  time,  the  problem  is 
much  more  easily  handled  by  the  new  method  and  the  oscillations 
about  the  exact  solution  are  much  reduced.  A  single  discontinuity 
in  velocity  has  been  replaced  by  two  discontinuities  in  acceleration. 
The  exact  solution  for  this  case  can  be  obtained  as  described  in 
Appendix  A.  The  computer  program  for  the  new  boundary  condition 
needs  only  to  be  modified  by  specifying  the  actual  motion. 

A  computation  was  made  for  the  case  of  a  finite  velocity-rise¬ 
time,  with  the  velocity  of  the  impacted  end  of  the  string  rising 
linearly  from  0  to  v0,  in  the  time  interval  0  to  .5*  This  rise  time 
was  typical  of  what  was  observed  in  much  of  the  experimental  data. 

In  Fig.  23  the  exact  and  computed  solutions  for  row  105  are  compared. 
The  usual  square  lattice  was  taken,  with  n  =  60j  hence  t  =  1.75* 

The  leading  edge  of  the  strain  wave  has  travelled  1.75  times  the 
length  of  the  string  and  the  trailing  edge  has  traveled  1.25  times 
the  length  of  the  string.  As  before,  the  straight  lines  represent 
the  exact  solution  and  the  curved  line  represents  the  computer 
solution.  Comparing  Fig.  23  with  Fig.  22  we  see  that  the  substit¬ 
ution  of  a  continuous  rise  in  velocity  for  the  instantaneous  rise 
has  very  greatly  reduced  the  tendency  of  the  computer  solution  to 
oscillate.  Whereas  in  Fig.  22  the  greatest  departure  of  the 
computer  strain  from  the  exact  solution  was  28  percent,  the  greatest 
departure  in  Fig.  23  is  only  3  percent. 

In  an  acutal  situation  the  finite  rise-time  of  velocity  will 
help  to  reduce  the  errors,  due  to  oscillations,  of  the  successive- 
substitution  computer- solution.  Delayed  response  of  the  material 
will  also  make  it  easier  for  oscillations  in  the  computer  solution 
to  be  avoided.  Nonlinearity  of  the  material  will  generally,  but 
not  always,  have  a  salutary  effect.  All  in  all,  we  consider  the 
oscillations  of  the  computer  solution  to  be  troublesome  but  not 
serious  and  accept  the  validity  of  the  computer  solution  as 
established. 

One  might  expect  the  oscillations  of  the  computer  solution  to 
decrease  as  the  lattice  spacing  h  is  decreased.  A  test  was  made, 
putting  n  =  1/h  =  20,  kO,  60,  and  100,  and  solving  our  zero-rise¬ 
time  problem.  Increasing  n  greatly  improves  (increases)  the  steep¬ 
ness  of  the  computed  wave  front.  It  reduces  the  wavelength  of  the 
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Fig.  23  Strain  as  obtained  from  digital  computer  calculations,  for  n  =  60  with  j  =  105  (t  =  1.75) 
In  earlier  graphs  the  acceleration  of  the  impacted  end  of  the  string  has  occurred  in  zero 
time  whereas  here  the  acceleration  is  finite  and  constant  and  occurs  in  the  interval 


oscillations,  but  it  does  not  reduce  the  amplitude  of  the  oscillations 
as  much  as  might  be  hoped. 

Fortunately,  the  amplitude  of  the  oscillations  in  a  particular 
region  falls  as  the  wave  front  travels  away  from  that  region.  This 
is  shown  in  Fig.  2b.  The  lower  curve  (t  =  .5)  represents  the  50th 
row  of  the  calculated  lattice  points.  The  upper  curve  (t  =  1.5) 
represents  the  150th  row  of  the  lattice.  In  the  right-hand  half  of 
the  string  there  has  been  in  the  meantime  no  change  in  the  exact 
solution.  The  successive-substitution  process  has,  however,  modi¬ 
fied  and  smoothed  the  oscillations  in  this  region.  The  wavelength 
of  the  oscillations  at  the  later  time  is  only  about  half  of  that  at 
the  earlier  time.  The  amplitude  of  the  oscillations  in  the  mean¬ 
time  has  been  reduced  by  a  factor  of  5  or  4.  These  two  isolated 
samplings  of  the  oscillations  do  not  of  course  show  how  the  oscil¬ 
lations  behave  at  intermediate  times,  but  it  is  reasonable  to 
assume  that  they  gradually  diminish  in  size. 

A  smooth  curve  faired  through  the  computer  solution,  but  not 
following  the  oscillations,  would  agree  with  the  exact  solution 
everywhere  except  in  the  immediate  neighborhood  of  the  front  of 
the  square  wave.  As  the  oscillations  diminish  with  time,  the 
computer  solution  settles  down  to  agree  more  and  more  closely  with 
the  exact  solution. 

This  completes  the  description  of  the  steps  taken  in  devising 
the  integral  equation  method  as  a  means  for  solving  partial 
differential  equations,  in  devising  the  modified  method  of  succes¬ 
sive  substitutions  for  solving  the  integral  equations,  and  in 
testing  the  solution  obtained  against  the  exact  solution  in  the 
case  of  the  linear  string. 

This  description  will  be  continued  in  the  second  and  third 
reports  of  this  series* in  which  integral  equation  methods  will 
be  applied  to  the  non-linear  and  simultaneous  partial  differential 
equations  which  arise  in  the  treatment  of  strings  composed  of 
nonlinear  and  time-dependent  materials. 
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Fig.  2h  Strain  calculated  with  n  =  100,  at  n  =  50  (t  =  .5)  and  n  =  150  (t  =  1*5)#  Note  that  in 
the  time  interval  between  the  two  curves,  the  oscillations  in  the  computer  solution  have 
been  greatly  reduced  by  the  successive- substitution  process.  Infinite  initial  acceleration 
(square  wave). 


APPENDIX  A 


Calculation  of  the  Motion  of  the  Linear  String  When  the  Motion 

of  the  Moving  End  is  Specified  Arbitrarily. 


In  the  original  treatment  the  velocity  vQ  of  the  right  hand  end 
of  the  string  was  considered  to  be  constant.  This  abrupt  change  in 
the  end  velocity  when  t  =  0  gives  rise  to  the  corners  in  the  surface 
s(x, t),  as  shown  in  Figs.  9  and  10.  These  corners  in  turn  were  dif¬ 
ficult  to  approximate  accurately  by  the  numerical  process,  and  for  a 
given  accuracy  would  require  the  mesh  of  the  lattice  used  to  be 
unduly  small.  Since  in  an  actual  experiment  the  end  velocity  cannot 
be  changed  abruptly  it  appears  that  the  above  difficulty  is  artifi¬ 
cial,  and  can  probably  be  eliminated  if  a  realistic  end-velocity- 
versus-time  function  is  used  instead  of  a  step  function.  If  this 
is  done  the  numerical  procedure  is  altered  only  in  that  sn .  is 
taken  from  the  given  motion  of  the  end,  and  is  no  longer  linear  in 
j.  The  modification  in  the  numerical  process  is  hence  trivial.  To 
test  the  accuracy  of  the  numerical  solution  obtained  under  this  new 
condition  it  is,  as  before,  desirable  to  consider  a  case  where  the 
exact  solution  can  be  obtained.  The  linear  string  provides  such  a 
case,  as  we  shall  now  see. 


Referring  to  (51),  which  pertains  to  the  linear  string,  we  see 
that  if  Si(x,t)  and  s2(x,t)  are  the  solutions  obtained,  respectively, 
for  the  given  initial  and  end  conditions 


s(x,o)  =  Si(x,o),  ||  = 


s(o,t)  =  Si(0,t), 
and  the  conditions 


t=o  t=o 
s(L,t)  =  Si(L,t); 


(1^2) 


s(x,o)  =  s2(x,o), 
s(o,t)  =  S2(o,t), 


5s  _  5s2 
dtt=0  5t  t=o 

s(L,t)  =  s2  (L,t); 


(143) 


then  the  integral  equation  (51)  is  satisfied  if  a  subscript  1  be 
affixed  to  s  wherever  it  appears.  Similarly  the  integral  equation 
is  satisfied  if  a  subscript  2  be  so  affixed.  Adding  the  two 
equations  so  obtained  we  get 


Si  +  s2  =  Si  (x,o)  +  s2(x,o)  +  t 


t=o 


a2  (Si  +  S2) 

— 


dt  dt. 


(144) 


which  equation  states  that  if  the  two  sets  of  initial  and  boundary 
conditions  (142)  and  (143)  be  superimposed,  then  the  solution  is 
obtained  by  superimposing  the  two  component  solutions  Si(x, t)  and 
s2(x,  t).  Superposition  is  therefore  permissible. 

In  those  cases  which  are  of  interest  to  us  s(o,t)  and  s(x,o) 
both  vanish,  and  (ds/^t)^._0  vanishes  if  0  <  x  <  L.  Our  problem  is 
to  find  s(x,  t)  when  the  initial  and  end  conditions  are  given  by 


s(x,o)  =  0, 
s(o,t)  =  0, 


(ft)'  ^ 


0  <  x  <  L, 
0  <  t, 

0  <  x  <  L, 

x  =  L. 


(145) 


Using  superposition  let  us  denote  by  su  (x, t)  the  solution 
corresponding  to  (145)  with  v0(t)  =  1  —  that  is,  a  step  function. 
The  desired  function  vQ(t)  may  be  regarded  to  be  composed  of  an 
infinite  number  of  step  functions  as  indicated  in  Fig.  25*  These 
start  at  different  times,  and  our  objective  is  to  find  their  cumu¬ 
lative  effect  at  time  t.  The  elementary  step  function  of  magnitude 
dvQ  which  begins  at  time  r  acts  during  a  time  (t-r),  and  hence  con¬ 
tributes  an  amount  s^x,  (t-r)]  dvQ  to  s(x,t).  Since  s(x,t)  is 
composed  of  the  sum  of  such  contributions  it  follows  that 


s(x,t) 


(t) 

su(x, (t-r)]  dvQ  . 


(146) 


% 

* 


\ 

\ 
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Fig.  25  Representation  of  a  continuous  velocity  function  vQ(t)  by  a 
series  of  step  functions. 


Changing  the  integration  variable  from  vQ  to  t,  this  becomes 


i* 


s(x,t)  Sufx, (t-r)]  V0  (r)  dT, 


(147) 


where  the  prime  indicates  differentiation.  Integrating 
(147)  becomes 


s(x,t)  = 


8u[x,(t-r)]  v0(r)J  v0(t) 


by  parts 
dr  (148) 


where  in  the  integrand  t  is  replaced  by  (t-r)  after  the  partial  dif¬ 
ferentiation  has  been  carried  out.  Since 


su(x,o)  =  0,  and  vo(0)  =  0 
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(149) 


the  first  term  on  the  right  hand  side  of  (148)  vanishes,  leaving 


s(x,t)  = 


T  /dsu(x,t)\ 

J0  V,dt  k=  t-T 


v0(t)  dT. 


(150) 


Here  x  is  a  constant  parameter  and  t  is  held  constant  during  the 
integration. 


Referring  to  (88),  (90),  and  Fig.  9,  we  see  that 


dsu(x,t) 

— 5t — 


if  -  (nL-x)  <  t  <  -  (nL+x) 


and  n  *  1,  3,  5,  7,  ; 


8su(x,t) 

“St - 


=  0  for  all  other  values  of  t. 


(151) 


Plotted  as  a  function  of  t  for  fixed  x,  dsu(x,t)/8t  appears  graph¬ 
ically  as  shown  in  Fig.  26. 


Fig.  26  Time  derivative  of  the  displacement  due  to  an  end  velocity 
which  is  a  unit  step  function. 


Laying  this  off  backward  from  time  t  we  obtain  a  graph  of 

[dsu(x,  t  )/3t]t=^._T.  This  and  the  other  factor  in  the  integrand  of  (150), 

namely,  vq(t),  are  shown  in  Fig,  27-  Since  s(x,t)  is  the  integral  of 


Fig.  27  Graph  showing  the  two  factors  in  the  integrand  of  s,  when 
the  displacement  is  given  by  a  superposition  integral. 

the  product  of  the  two  factors  we  that  s(x,t)  is  given  by  the  total 
area  under  those  parts  of  the  curve  of  v0(t)  where  the  graph  of  the 
other  factor  [dsu(x,t)/dt]^._t_T  does  not  vanish.  As  t  increases, 

this  latter  graph  moves  to  the  right,  whereas  the  curve  of  v0(r) 
remains  stationary.  Since 


v0(T )  dr  =  s(L,Ta)-s(L,Tx)*&s(L,t) 


(152) 


we  see  that  s(x, t)  is  the  sum  of  the  increments  of  the  end  displace¬ 
ment  curve  s(L,r)  which  occur  over  the  intervals  in  which 
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Fig.  28  Graphs  of  functions  to  be  used  in  finding  s(x, t). 
&su(x,t)/dt]  does  not  vanish,  as  indicated  in  Fig.  28. 

Z-Z-T 

end  displacement  curve  s(L,r)  is,  of  course,  given. 


The 


APPENDIX  B 


Calculation  of  the  Numerical  Approximation  of  the  Motion  of 

the  Linear  String  When  the  Motion  of  the  Moving  End  is 

Specified  Arbitrarily. 

if 

Let  us  suppose  that  the  displacements  s^j  at  the  various  lattice 
points  involved  in  the  approximate  solution  for  the  motion  of  the 
linear  string  when  the  moving  end  has  a  constant  unit  velocity  are 
given.  We  wish  to  determine  the  displacements  s^j  when  the  moving 
end  has  the  arbitrarily  specified  motion  snj.  Placing 

vnj  “*h^snj  “  n(sn,j+i"  8n j ^  J  =  ^  ? >  •••  (!53) 

i.1. 

we  see  that  the  velocity  in  the  (j+l)  time  subinterval  at  the 
moving  end,  which  is  the  (j+l)th  subinterval  in  the  right  hand  col¬ 
umn  of  lattice  points  is  vQj,  Fig.  29.  As  indicated  in  the  figure. 


Fig.  29  Velocity  of  moving  end  of  string.  The  actual  velocity  is 

approximated  by  adding  a  finite  increment  to  vn  at  each 
successive  row. 


the  velocity  of  the  moving  end  may  be  regarded  to  be  composed  of  a 
number  of  abrupt  increments  which  occur  at  different  times,  but  con¬ 
tinue  indefinitely  after  they  do  occur.  Placing 


Av, 


nJ 


=  v, 


n,J+i  "  n,j 


(15*0 


we  see  that  the  increase  in  velocity  at  the  time  corresponding  to  the 

value  k  of  the  index  is  Av  , 

n,  k-i 

Since  the  system  is  linear  the  displacement,  s-^j  may  be  regarded 
to  be  composed  of  the  effects  of  the  various  velocity  increments, 
these  being  superimposed.  The  contribution  to  s^  of  the  velocity 
increment  Av  .  which  occurs  at  k  is 


si, j-k  Avn,k-i  ; 


hence  adding  these  contributions  we  obtain  finally 


(155) 


Starting  with  the  given  values  of  snj  the  Av's  in  (156)  can  easily  be 
obtained  by  making  a  table  of  first  and  second  differences. 

If  s^j  is  obtained  from  the  solution  obtained  by  solving  dif¬ 
ference  equations,  then  (156)  gives  the  corresponding  solution  when 
the  motion  of  the  moving  end  is  specified  arbitrarily.  This  can  be 
used  to  check  the  result  obtained  by  the  successive  approximation 
process,  which  process  provides  for  the  arbitrary  specification  of 
the  motion  of  the  moving  end. 

If  the  specified  motion  of  the  moving  end  has  no  discontinuity 
in  velocity,  the  approximate  solution  s^j  may  be  very  much  more 
accurate  than  .,  which  involves  such  a  discontinuity.  This  fact 
does  not  affect  £he  validity  of  (156),  which  gives  s^  using  s*^. 
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